GEANT321/gcons/gplmat.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:15  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/04 10/01/95  17.36.36  by  S.Ravndal
  11 *-- Author :
  12       SUBROUTINE GPLMAT(IMATES,IPART,MECAN,KDIN,TKIN,IDM)
  13 C.
  14 C     ******************************************************************
  15 C.    *                                                                *
  16 C.    *       INTERPOLATE and PLOT  the DE/DX and Cross sections       *
  17 C.    *       tabulated in JMATE banks corresponding to :              *
  18 C.    *       material IMATE, particle IPART, mecanism name HMECAN,    *
  19 C.    *       kinetic energies TKIN                                    *
  20 C.    *                                                                *
  21 C.    *      The MECANism name can be :                                *
  22 C.    *      'HADF'   'INEF'   'ELAF'   'FISF'   'CAPF'                *
  23 C.    *      'HADG'   'INEG'   'ELAG'   'FISG'   'CAPG'                *
  24 C.    *      'LOSS'   'PHOT'   'ANNI'   'COMP'   'BREM'                *
  25 C.    *      'PAIR'   'DRAY'   'PFIS'   'RAYL'   'HADG'                *
  26 C.    *      'MUNU'   'RANG'   'STEP'                                  *
  27 C.    *                                                                *
  28 C.    *       For Hadronic particles it also computes the              *
  29 C.    *       hadronic cross section from FLUKA ( '***F' ) or          *
  30 C.    *       GHEISHA ( '***G' ) programs:                             *
  31 C.    *       HADF or HADG -- total                                    *
  32 C.    *       INEF or INEG -- inelastic                                *
  33 C.    *       ELAF or ELAG -- elastic                                  *
  34 C.    *       FISF or FISG -- fission (0.0 for FLUKA)                  *
  35 C.    *       CAPF or CAPG -- neutron capture (0.0 for FLUKA)          *
  36 C.    *                                                                *
  37 C.    *             Input parameters                                   *
  38 C.    *  IMATE   Geant material number                                 *
  39 C.    *  IPART   Geant particle number                                 *
  40 C.    *  MECAN  mechanism name of the bank to be fetched               *
  41 C.    *  KDIM   dimension of the arrays TKIN , VALUE                   *
  42 C.    *  TKIN   array of kinetic energy of incident particle (in Gev)  *
  43 C.    *  IDM    convention for histogramming mode :                    *
  44 C.    *         IDM.gt.0  fill , print , keep histogram(s)             *
  45 C.    *         IDM.eq.0  fill , print , delete histogram(s)           *
  46 C.    *         IDM.lt.0  fill , noprint , keep histogram(s)           *
  47 C.    *           The histogram IDentificator will be :                *
  48 C.    *             10000*imate + 100*ipart + imeca                    *
  49 C.    *          where IMECA is the link number in stucture JMATE      *
  50 C.    *          (see Geant3 writeup CONS 199)                         *
  51 C.    *           for 'HADG'  imeca = 17                               *
  52 C.    *                                                                *
  53 C.    *    ==>Called by : <USER>                                       *
  54 C.    *       Authors    R.Brun, M.Maire    *********                  *
  55 C.    *                                                                *
  56 C.    ******************************************************************
  57 C.
  58 #include "geant321/gcbank.inc"
  59 #include "geant321/gcnum.inc"
  60 #include "geant321/gconsp.inc"
  61 #include "geant321/gcunit.inc"
  62 #include "geant321/gcphys.inc"
  63       PARAMETER (MMX= 201,NCOL= 5)
  64       CHARACTER*(*) MECAN
  65       CHARACTER*4 MECA , KU(NCOL)
  66       CHARACTER   NAPART*16 , NAMATE*16 , CHTITL*68
  67       DIMENSION   TKIN(KDIN),VALUE(MMX),SIGT(MMX),PCUT(5)
  68       DIMENSION   KI(NCOL),EK(NCOL)
  69       LOGICAL     LXBARN
  70 *
  71 #include "geant321/gcnmec.inc"
  72 *
  73 *     ------------------------------------------------------------------
  74 *
  75       IF (KDIN.LE.0)  GO TO 999
  76       KDIM = MIN(KDIN,MMX)
  77       IF(IMATES.LT.0) THEN
  78          LXBARN=.TRUE.
  79       ELSEIF(IMATES.GT.0) THEN
  80          LXBARN=.FALSE.
  81       ELSE
  82          GOTO 999
  83       ENDIF
  84       IMATE=ABS(IMATES)
  85 *
  86       IF (JMATE.LE.0) GO TO 999
  87       IF (IMATE.GT.NMATE) GO TO 80
  88       JMA = LQ(JMATE-IMATE)
  89       IF  (JMA.LE.0) GO TO 80
  90       CALL UHTOC(IQ(JMA+1),4,NAMATE,16)
  91       IF(LXBARN) THEN
  92          CMIBAR=Q(JMA+6)/(AVO*Q(JMA+8))
  93       ELSE
  94          CMIBAR=1.
  95       ENDIF
  96 *
  97       IF (JPART.LE.0) GO TO 999
  98       IF (IPART.LE.0) GO TO 999
  99       IF (IPART.GT.NPART) GO TO 80
 100       JPA = LQ(JPART-IPART)
 101       IF (JPA.LE.0) GO TO 80
 102       CALL UHTOC(IQ(JPA+1),4,NAPART,16)
 103 *
 104 * *** Print  bin meaning
 105       IF (IDM.GE.0) THEN
 106          CHMAIL='1'
 107          CALL GMAIL(0,0)
 108          CHMAIL=' '
 109          CHMAIL(31:)='Kinetic energy bin meaning'
 110          CALL GMAIL(0,0)
 111          CHMAIL(31:)='--------------------------'
 112          CALL GMAIL(0,1)
 113          NROW = (KDIM-1)/NCOL + 1
 114          DO 20  IR=1,NROW
 115             DO 10  IC=1,NCOL
 116                IKB = IR + (IC-1)*NROW
 117                IF (IKB.GT.KDIM) IKB=KDIM
 118                KI(IC) = IKB
 119                CALL GEVKEV(TKIN(IKB),EK(IC),KU(IC))
 120    10       CONTINUE
 121             WRITE(CHMAIL,10200) (KI(IC),EK(IC),KU(IC),IC=1,NCOL)
 122             CALL GMAIL(0,0)
 123    20    CONTINUE
 124       ENDIF
 125 *
 126       BIGINV= 1000./BIG
 127       DO 30  JMX = 1, MMX
 128          SIGT(JMX) = 0.
 129    30 CONTINUE
 130       IF(MECAN.EQ.'ALL') THEN
 131          N1 = 1
 132          N2 = NMECA
 133       ELSE
 134          N1 = 0
 135          DO 40  IMECA=1,NMECA
 136             IF(MECAN.EQ.CHNMEC(IMECA)) THEN
 137                N1 = IMECA
 138             ENDIF
 139    40    CONTINUE
 140          IF(N1.EQ.0) THEN
 141             WRITE(CHMAIL,'('' *** GPLMAT: Mechanism '',A,
 142      +      '' not implemented'')') MECAN
 143             CALL GMAIL(0,0)
 144             GOTO 999
 145          ENDIF
 146          N2 = N1
 147       ENDIF
 148       DO 60  IMEC = N1,N2
 149 C
 150          IF (MECAN.EQ.'ALL') THEN
 151              IF (CHNMEC(IMEC).EQ.'RANG') GO TO 60
 152              IF (CHNMEC(IMEC).EQ.'STEP') GO TO 60
 153          END IF
 154 C
 155          IF(CHNMEC(IMEC).NE.'NULL') THEN
 156             MECA = CHNMEC(IMEC)
 157             CALL GFTMAT(IMATE,IPART,MECA,KDIM,TKIN,VALUE,PCUT,IXST)
 158             IF(IXST.EQ.0) GO TO 60
 159 *
 160 * ***    Book histogram
 161             ISIG = 0
 162             IF (MECA.EQ.'LOSS') THEN
 163                CHTITL = NAPART//' in '//NAMATE//'   dE/dx (MeV/cm)'
 164             ELSEIF (MECA.EQ.'RANG') THEN
 165                CHTITL = NAPART//' in '//NAMATE//'   Stopping range (cm)'
 166             ELSEIF (MECA.EQ.'STEP') THEN
 167                CHTITL = NAPART//' in '//NAMATE//'   continuous step '
 168      +         //'(cm)'
 169             ELSE
 170                CHTITL = NAPART//' in '//NAMATE//'   '//MECA// ' cross '
 171      +         //'section'
 172                IF(LXBARN) THEN
 173                   CHTITL(LNBLNK(CHTITL)+1:) = ' (barn)'
 174                ELSE
 175                   CHTITL(LNBLNK(CHTITL)+1:) = ' (1/cm)'
 176                ENDIF
 177                ISIG = 1
 178             ENDIF
 179 *
 180             ID = 10000*IMATE + 100*IPART + IMEC
 181             CALL HBOOKB(ID,CHTITL,KDIM-1,TKIN,0.)
 182 *
 183 * ***    Fill histogram
 184 *
 185             VALMI = MAX (BIGINV,VMAX(VALUE,KDIM)*1.E-8)
 186             DO 50  IKB = 1,KDIM
 187                IF (MECA.NE.'LOSS'.AND.MECA.NE.'RANG'
 188      +             .AND.MECA.NE.'STEP')
 189      +             VALUE(IKB)=VALUE(IKB)*CMIBAR
 190                IF (VALUE(IKB).GE.VALMI) THEN
 191                   CALL HFILL(ID,TKIN(IKB),0.,VALUE(IKB))
 192                ENDIF
 193                IF (ISIG.EQ.1) THEN
 194                   IF(MECA(1:3).NE.'INE'.AND.MECA(1:3).NE.'ELA'.AND.
 195      +            MECA(1:3).NE.'FIS'.AND.MECA(1:3).NE.'CAP'.AND.
 196      +            MECA(1:3).NE.'HAD'.AND.IMEC.LT.IBLOWN) THEN
 197                      SIGT(IKB) = SIGT(IKB) + VALUE(IKB)
 198                   ELSE IF (MECA(1:3).EQ.'HAD') THEN
 199                      IF ((MECA.EQ.'HADG'.AND.IHADR.LE.2).OR. (MECA.EQ.
 200      +               'HADF'.AND.IHADR.EQ.4)) THEN
 201                         SIGT(IKB) = SIGT(IKB) + VALUE(IKB)
 202                      ENDIF
 203                   ENDIF
 204                ENDIF
 205    50       CONTINUE
 206             CALL HIDOPT(ID,'LOGY')
 207             IF(IDM.GE.0) CALL HPHIST(ID,' ',0)
 208             IF(IDM.EQ.0) CALL HDELET(ID)
 209          ENDIF
 210    60 CONTINUE
 211 *
 212 * *** plot total cross section and mean free path
 213       IF (MECAN.EQ.'ALL') THEN
 214          CHTITL= NAPART//' in '//NAMATE//'   total cross section'
 215          IF(LXBARN) THEN
 216             CHTITL(LNBLNK(CHTITL)+1:) = ' (barn)'
 217          ELSE
 218             CHTITL(LNBLNK(CHTITL)+1:) = ' (1/cm)'
 219          ENDIF
 220          ID = 10000*IMATE + 100*IPART + NMECA+1
 221          CALL HBOOKB(ID,CHTITL,KDIM-1,TKIN,0.)
 222 *
 223          CHTITL= NAPART//' in '//NAMATE//'   total mean free path (cm)'
 224          II = ID + 1
 225          CALL HBOOKB(II,CHTITL,KDIM-1,TKIN,0.)
 226 *
 227          VALMI = MAX (BIGINV,VMAX( SIGT,KDIM)*1.E-8)
 228          DO 70  IKB = 1,KDIM
 229             IF (SIGT(IKB).GE.VALMI) THEN
 230                CALL HFILL(ID,TKIN(IKB),0.,       SIGT(IKB))
 231                CALL HFILL(II,TKIN(IKB),0.,CMIBAR/SIGT(IKB))
 232             ENDIF
 233    70    CONTINUE
 234          CALL HIDOPT(ID,'LOGY')
 235          IF(IDM.GE.0) CALL HPHIST(ID,' ',0)
 236          IF(IDM.EQ.0) CALL HDELET(ID)
 237 *
 238          CALL HIDOPT(II,'LOGY')
 239          IF(IDM.GE.0) CALL HPHIST(II,' ',0)
 240          IF(IDM.EQ.0) CALL HDELET(II)
 241       ENDIF
 242 *
 243       GO TO 999
 244 *
 245    80 WRITE(CHMAIL,10000) IMATE ,IPART
 246       CALL GMAIL(0,0)
 247 10000 FORMAT(' ***** GPLMAT error : material',I4,
 248      +       '  or particle',I4,' not defined'   )
 249 10100 FORMAT(6X,'BCUTE =',F6.2,A4,3X,'BCUTM =',F6.2,A4,3X,
 250      +             'DCUTE =',F6.2,A4,3X,'DCUTM =',F6.2,A4,3X,
 251      +            'PPCUTM =',F6.2,A4 )
 252 10200 FORMAT(1X,5('   bin ',I3,' =',F7.2,A4))
 253   999 END