GEANT321/gphys/gphxsi.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:31  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.22  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GPHXSI
  13 *.
  14 *.    ******************************************************************
  15 *.    *  Creates PHXS bank containing x-section constants              *
  16 *.    *                                                                *
  17 *.    *    ==>CALLED BY : GPHYSI                                       *
  18 *.    *       AUTHOR : J. CHWASTOWSKI                                  *
  19 *.    *                                                                *
  20 *.    ******************************************************************
  21 *.
  22 *      The photoelectric effect x-section bank
  23 *       1                 - Number of elements neded to create medium = NZ
  24 *       2      <-> NZ+1   -  Z of elements
  25 *       NZ+1   <-> 2*NZ+1 -  0
  26 *       2+2*NZ <-> 3*NZ+1 -  weight of x-section constants
  27 *       3*NZ+2 <-> top    -  X-section constants
  28 #include "geant321/gcbank.inc"
  29 #include "geant321/gcphys.inc"
  30 #include "geant321/gcjloc.inc"
  31 #include "geant321/gconsp.inc"
  32 #include "geant321/gcunit.inc"
  33 #include "geant321/gcpmxz.inc"
  34 #include "geant321/gcphxs.inc"
  35 #include "geant321/gc10ev.inc"
  36       CHARACTER*1 CHSHEL, CHGROU
  37       DIMENSION EUP(MAXINT),TMP(MAXPOW,MAXINT),ESHL(24)
  38       PARAMETER (NXSB=73,NXSBF=70)
  39       Z = Q(JMA+7)
  40       IPHXSP = 0
  41 * Check Z range of validity for Sandia parametrization
  42       IF(Z.GE.1.AND.Z.LE.MAXELZ) THEN
  43 *
  44 * Find number of elements neded to create current medium
  45 c
  46 * Is this medium a mixture ?
  47          NMIX = Q(JMA+11)
  48          IF(NMIX.GT.1) THEN
  49             NZ = 0
  50             JMIXT=LQ(JMA-5)
  51             DO 10 I = 1,NMIX
  52                ZCUR = Q(JMIXT+NMIX+I)
  53                IF(ZCUR.NE.INT(ZCUR)) THEN
  54 *
  55 * When Z is non integer you need 2 elements
  56                   NZ=NZ+2
  57                ELSE
  58                   NZ=NZ+1
  59                ENDIF
  60    10       CONTINUE
  61             CALL GWORK(3*NZ+1)
  62             WS(1) = NZ
  63 *
  64 * Calculate weigths
  65 *
  66             K = 1
  67             DO 20 I = 1,NMIX
  68                ZCUR = Q(JMIXT+NMIX+I)
  69                HZCUR = INT(ZCUR)
  70                WS(1+K) = HZCUR
  71                IF(ZCUR.NE.HZCUR) THEN
  72                   WS(1+2*NZ+K) = (HZCUR+1.-ZCUR)*Q(JMIXT+2*NMIX+I)
  73                   K = K+1
  74                   WS(1+K) = HZCUR+1
  75                   WS(1+2*NZ+K) = (ZCUR-HZCUR)*Q(JMIXT+2*NMIX+I)
  76                ELSE
  77                   WS(1+2*NZ+K) = Q(JMIXT+2*NMIX+I)
  78                ENDIF
  79                K = K+1
  80    20       CONTINUE
  81 *
  82 * Do Z values repeat ?
  83 *
  84             K = NZ
  85             DO 40 I = 1,NZ-1
  86                Z1 = WS(1+I)
  87                IF(Z1.GT.0.0) THEN
  88                   DO 30 J = I+1,NZ
  89                      Z2 = WS(1+J)
  90                      IF(Z1.EQ.Z2) THEN
  91 *                        WS(1+NZ+I) = WS(1+NZ+I)+WS(1+NZ+J)
  92                         WS(1+2*NZ+I) = WS(1+2*NZ+I)+WS(1+2*NZ+J)
  93                         WS(1+J) = -WS(1+J)
  94                         K = K-1
  95                      ENDIF
  96    30             CONTINUE
  97                ENDIF
  98    40       CONTINUE
  99 * Now you can book a pht. eff. x-sec. constant bank and hang it as
 100 * a first struc. link from JPHOT.
 101 * From this bank you hang the banks for each separate Z !!!
 102             NW = NZ*NXSB+2
 103             IF(K.EQ.1) NW = NXSB+2
 104             CALL MZBOOK(IXCONS,JPHXS,JPHOT,-1,'PHXS',NZ,NZ,NW,3,0)
 105 * Fill JPHXS bank and calculates the weights
 106             Q(JPHXS+1) = K
 107             NZN = K
 108             K = 1
 109             DO 50 I = 1,NZ
 110                IF(WS(1+I).GT.0.0) THEN
 111                   Q(JPHXS+1+K) = WS(1+I)
 112 *                  Q(JPHXS+1+K+NZN) = WS(1+I+NZ)
 113                   Q(JPHXS+1+K+2*NZN) = WS(1+I+2*NZ)
 114                   K = K+1
 115                ENDIF
 116    50       CONTINUE
 117          ELSE
 118 * Current medium consists of one "element"
 119             NZ = 1
 120             ZCUR = Z
 121             HZCUR = INT(ZCUR)
 122             IF(MOD(ZCUR,HZCUR).EQ.0) THEN
 123 * Now you can book a pht. eff. x-sec. constant bank and hang it as
 124 * a first struc. link from JPHOT.
 125                NW = NXSB+2
 126                CALL MZBOOK(IXCONS,JPHXS,JPHOT,-1,'PHXS',1,1,NW,3,0)
 127                Q(JPHXS+1) = 1.
 128                Q(JPHXS+2) = ZCUR
 129                Q(JPHXS+4) = 1.
 130             ELSE
 131 * Somebody is cheating. We need two elements
 132                NZ = NZ+1
 133                NW = 2*NXSB+2
 134                CALL MZBOOK(IXCONS,JPHXS,JPHOT,-1,'PHXS',2,2,NW,3,0)
 135                Q(JPHXS+1) = 2.
 136                Q(JPHXS+2) = HZCUR
 137                Q(JPHXS+3) = HZCUR+1
 138                Q(JPHXS+6) = (HZCUR+1.-ZCUR)
 139                Q(JPHXS+7) = (ZCUR-HZCUR)
 140             ENDIF
 141          ENDIF
 142 * We passed the bank booking phase and we can start real work
 143          NZ = Q(JPHXS+1)
 144 * Create a temporary working space
 145          IBINS = 5*NZ+NZ*(MAXPOW+1)*(MAXINT+1)
 146          CALL GWORK(IBINS)
 147          DO 60 I = 1,IBINS
 148             WS(I) = 0.0
 149    60    CONTINUE
 150          KS = 5*NZ
 151          L = 0
 152          DO 170 N = 1,NZ
 153             JPHXS = LQ(JPHOT-1)
 154             ZCUR = Q(JPHXS+1+N)
 155             IZ = ZCUR
 156             DO 70 I = 1,24
 157                ESHL(I) = 0.0
 158    70       CONTINUE
 159             CALL GFSHLS(ZCUR,ESHL,NSHL)
 160             DO 80 I = 1,NSHL
 161                ESHL(I) = ESHL(I)*1.E-3
 162    80       CONTINUE
 163 * Use Sandia data.
 164 * Find out the interval upper limit.
 165             IMAX = 0
 166             DO 90 I = 1,MAXINT
 167                CHSHEL=CRNGUP(I,IZ)(1:1)
 168                CHGROU=CRNGUP(I,IZ)(2:2)
 169                IF(CHSHEL.EQ.'I') THEN
 170                   EUP(I) = 55.E11
 171                   IMAX = I
 172                   GO TO 100
 173                ELSEIF(CHSHEL.EQ.'K') THEN
 174                   EUP(I) = ESHL(1)
 175                ELSEIF(CHSHEL.EQ.'L') THEN
 176                   IF(CHGROU.EQ.'1') THEN
 177                      EUP(I) = ESHL(2)
 178                   ELSEIF(CHGROU.EQ.'2') THEN
 179                      EUP(I) = ESHL(3)
 180                   ELSEIF(CHGROU.EQ.'3') THEN
 181                      EUP(I) = ESHL(4)
 182                   ELSE
 183                      CHMAIL = ' GPHXSI Error: Inconsistent data for L '
 184      +               //'shell.  Maximum coded L3. Found shell '
 185      +               //'name: '//CRNGUP(I,IZ)
 186                      CALL GMAIL(0,0)
 187                   ENDIF
 188                ELSEIF(CHSHEL.EQ.'M') THEN
 189                   IF(CHGROU.EQ.'1') THEN
 190                      EUP(I) = ESHL(5)
 191                   ELSEIF(CHGROU.EQ.'2') THEN
 192                      EUP(I) = ESHL(6)
 193                   ELSEIF(CHGROU.EQ.'3') THEN
 194                      EUP(I) = ESHL(7)
 195                   ELSEIF(CHGROU.EQ.'4') THEN
 196                      EUP(I) = ESHL(8)
 197                   ELSEIF(CHGROU.EQ.'5') THEN
 198                      EUP(I) = ESHL(9)
 199                   ELSE
 200                      CHMAIL = ' GPHXSI Error: Inconsistent data for M '
 201      +               //'shell.  Maximum coded M5. Found shell '
 202      +               //'name: '//CRNGUP(I,IZ)
 203                      CALL GMAIL(0,0)
 204                   ENDIF
 205                ELSEIF(CHSHEL.EQ.'N') THEN
 206                   IF(CHGROU.EQ.'1') THEN
 207                      EUP(I) = ESHL(10)
 208                   ELSEIF(CHGROU.EQ.'2') THEN
 209                      EUP(I) = ESHL(11)
 210                   ELSEIF(CHGROU.EQ.'3') THEN
 211                      EUP(I) = ESHL(12)
 212                   ELSEIF(CHGROU.EQ.'4') THEN
 213                      EUP(I) = ESHL(13)
 214                   ELSEIF(CHGROU.EQ.'5') THEN
 215                      EUP(I) = ESHL(14)
 216                   ELSEIF(CHGROU.EQ.'6') THEN
 217                      EUP(I) = ESHL(15)
 218                   ELSEIF(CHGROU.EQ.'7') THEN
 219                      EUP(I) = ESHL(16)
 220                   ELSE
 221                      CHMAIL = ' GPHXSI Error: Inconsistent data for N '
 222      +               //'shell. Maximum coded N7. Found shell name:'
 223      +               //' '//CRNGUP(I,IZ)
 224                      CALL GMAIL(0,0)
 225                   ENDIF
 226                ELSEIF(CHSHEL.EQ.'O') THEN
 227                   IF(CHGROU.EQ.'1') THEN
 228                      EUP(I) = ESHL(17)
 229                   ELSEIF(CHGROU.EQ.'2') THEN
 230                      EUP(I) = ESHL(18)
 231                   ELSEIF(CHGROU.EQ.'3') THEN
 232                      EUP(I) = ESHL(19)
 233                   ELSEIF(CHGROU.EQ.'4') THEN
 234                      EUP(I) = ESHL(20)
 235                   ELSEIF(CHGROU.EQ.'5') THEN
 236                      EUP(I) = ESHL(21)
 237                   ELSE
 238                      CHMAIL = ' GPHXSI Error: Inconsistent data for O '
 239      +               //'shell.  Maximum code O5. Found shell name:'
 240      +               //' '//CRNGUP(I,IZ)
 241                      CALL GMAIL(0,0)
 242                   ENDIF
 243                ELSEIF(CHSHEL.EQ.'P') THEN
 244                   IF(CHGROU.EQ.'1') THEN
 245                      EUP(I) = ESHL(22)
 246                   ELSEIF(CHGROU.EQ.'2') THEN
 247                      EUP(I) = ESHL(23)
 248                   ELSEIF(CHGROU.EQ.'3') THEN
 249                      EUP(I) = ESHL(24)
 250                   ELSE
 251                      CHMAIL = ' GPHXSI Error: Inconsistent data for P '
 252      +               //'shell.  Maximum coded P3. Found shell '
 253      +               //'name: '//CRNGUP(I,IZ)
 254                      CALL GMAIL(0,0)
 255                   ENDIF
 256                ELSE
 257                   READ(CRNGUP(I,IZ),'(F6.0)') EUP(I)
 258                ENDIF
 259                IF(EUP(I).EQ.0.0) THEN
 260                   WRITE(CHMAIL,'(A44,I5)') ' GPHXSI Error: Upper limit '
 261      +            //'= 0. Interval:',I
 262                   CALL GMAIL(0,0)
 263                   STOP 14
 264                ENDIF
 265    90       CONTINUE
 266   100       CONTINUE
 267             DO 120 I = 1,IMAX
 268                DO 110 J = 1,MAXPOW
 269                   TMP(J,I) = COFS(J,I,IZ)*Q(JPHXS+1+2*NZ+N)
 270   110          CONTINUE
 271   120       CONTINUE
 272 * Copy upper limits and coofficients to a work bank
 273             K = KS+1
 274             WS(K) = GPOMIN(IZ)
 275             WS(K+1) = 0.0
 276             WS(K+2) = 0.0
 277             WS(K+3) = 0.0
 278             WS(K+4) = 0.0
 279             K = K+5
 280             DO 140 I = 1,IMAX
 281                WS(K) = EUP(I)
 282                K = K+1
 283                DO 130 J = 1,MAXPOW
 284                   WS(K) = TMP(J,I)
 285                   K = K+1
 286   130          CONTINUE
 287   140       CONTINUE
 288             KS = KS+(MAXINT+1)*(MAXPOW+1)
 289             IWS(NZ+N) = IMAX
 290             IWS(2*NZ+N) = 0
 291             IWS(3*NZ+N) = 0
 292 *            WS(N) = Q(JPHXS+1+2*NZ+N)
 293 * Create element x-section & final state bank
 294             CALL MZBOOK(IXCONS,JPHFN,JPHXS,-N,'PHFN',0,0,IMAX*5+1,3,0)
 295             Q(JPHFN+1) = IMAX
 296 * Update pointer
 297             JPHFN = JPHFN+1
 298             KFN = 1
 299 * Copy energy & x-section parameters
 300             DO 160 I = 1,IMAX
 301                Q(JPHFN+KFN) = EUP(I)
 302                KFN = KFN+1
 303                DO 150 J = 1,MAXPOW
 304 c                  Q(JPHFN+KFN) = TMP(J,I)*Q(JPHXS+1+2*NZ+N)
 305                   Q(JPHFN+KFN) = TMP(J,I)
 306                   KFN = KFN+1
 307   150          CONTINUE
 308   160       CONTINUE
 309 * Get shells decay parameters
 310             CALL GFSHDC(N,ZCUR)
 311   170    CONTINUE
 312 * Now find intervals and calculate the coofs for each
 313          K = 0
 314          JPHXS = LQ(JPHOT-1)
 315          IF(NZ.LT.2) THEN
 316 * Simple element so life is easy and nice
 317             Q(JPHXS+5) = IMAX+1
 318             JPHXS6 = JPHXS+6
 319             Q(JPHXS6) = GPOMIN(IZ)
 320             Q(JPHXS6+1) = 0.0
 321             Q(JPHXS6+2) = 0.0
 322             Q(JPHXS6+3) = 0.0
 323             Q(JPHXS6+4) = 0.0
 324             JPHXS6 = JPHXS6+5
 325             DO 190 I = 1,IMAX
 326                Q(JPHXS6+K) = EUP(I)
 327                K = K+1
 328                DO 180 J = 1,MAXPOW
 329                   Q(JPHXS6+K) = TMP(J,I)
 330                   K = K+1
 331   180          CONTINUE
 332   190       CONTINUE
 333             NPSH = NW-5*(IMAX+1)-5
 334          ELSE
 335 * More elements. It will not be so easy
 336 * The following code is difficult and probably there are better solutions
 337 * but I could think only about this one.
 338             IPHXSP = JPHXS+2+3*NZ
 339             IPIMAX = JPHXS+2+3*NZ
 340             IOFFST = (MAXPOW+1)*(MAXINT+1)
 341             DO 250 II = 1,NZ*(MAXINT+1)
 342 * Find the interval for which the upper limit is the smallest
 343                AMINV = 1.E20
 344                DO 200 I = 1,NZ
 345                   K = IWS(2*NZ+I)
 346                   IPOINT = 5*NZ+1+(I-1)*IOFFST+K
 347                   IEUP = 4*NZ+I
 348                   WS(IEUP) = WS(IPOINT)
 349                   IF(WS(IEUP).LT.AMINV) AMINV = WS(IEUP)
 350   200          CONTINUE
 351                L = 0
 352                DO 210 I = 1,NZ
 353                   IF(WS(4*NZ+I).LE.AMINV) L = L+1
 354   210          CONTINUE
 355                IF(L.LT.1) THEN
 356                   CHMAIL = ' GPHXSI Error: Zero intervals found.'
 357                   CALL GMAIL(0,0)
 358                   STOP 16
 359                ENDIF
 360 * Copy to JPHXS bank
 361                Q(IPIMAX) = Q(IPIMAX)+1.
 362                Q(IPHXSP+1) = AMINV
 363                IPHXSP = IPHXSP+1
 364                DO 230 J = 1,MAXPOW
 365                   QS = 0.0
 366                   DO 220 I = 1,NZ
 367                      K = IWS(2*NZ+I)
 368                      IPOINT = 5*NZ+1+(I-1)*IOFFST+K+J
 369 c                     QS = QS+WS(I)*WS(IPOINT)
 370                      QS = QS+WS(IPOINT)
 371   220             CONTINUE
 372                   Q(IPHXSP+1) = QS
 373                   IPHXSP = IPHXSP+1
 374   230          CONTINUE
 375                IF(L.EQ.NZ.AND.AMINV.EQ.55.E11) GO TO 260
 376 * Update local pointers
 377                DO 240 I = 1,NZ
 378                   IEUP = 4*NZ+I
 379                   IF(WS(IEUP).LE.AMINV) THEN
 380                      INZ2 = 2*NZ+I
 381                      IF(IWS(3*NZ+I).LT.IWS(I+NZ)) THEN
 382                         IWS(INZ2) = IWS(INZ2)+5
 383                         IWS(3*NZ+I) = IWS(3*NZ+I)+1
 384                      ENDIF
 385                   ENDIF
 386   240          CONTINUE
 387   250       CONTINUE
 388 * It is THE END of th x-secs. part, however you may not believe it.
 389   260       CONTINUE
 390             NIT=Q(IPIMAX)
 391             IIT = -1
 392   261       IIT = IIT+1
 393 *
 394 * It just may happen that some of the energy limits are the same
 395                IENE=IIT*5+1
 396                IF(ABS(Q(IPIMAX+IENE)-Q(IPIMAX+IENE+5)).LT.
 397      +                                Q(IPIMAX+IENE)*5E-4) THEN
 398                   DO 262 II=1,(NIT-IIT-1)*5
 399                      Q(IPIMAX+IENE+II)=Q(IPIMAX+IENE+II+5)
 400   262             CONTINUE
 401                   NIT=NIT-1
 402                   IIT=IIT-1
 403                ENDIF
 404             IF(IIT.LT.NIT-1) GO TO 261
 405             Q(IPIMAX)=NIT
 406             NPSH = NW-Q(IPIMAX)*5-3*NZ-2
 407          ENDIF
 408 * Return unused locations in JPHXSI bank to the system
 409          IF(NPSH.GT.0) THEN
 410             JPHXS = LQ(JPHOT-1)
 411             CALL MZPUSH(IXCONS,JPHXS,0,-NPSH,'R')
 412          ENDIF
 413       ELSEIF(Z.GT.100.) THEN
 414 * Just in case we got called
 415          CHMAIL = ' GPHXSI Error: Z > 100. No Sandia parameters. '
 416          CALL GMAIL(0,0)
 417       ENDIF
 418       END