5 * Revision 1.3 1996/09/30 14:26:05 ravndal
6 * Windows NT related modifications
8 * Revision 1.2 1996/04/26 12:21:12 ravndal
9 * NAP character*20 declaration included
11 * Revision 1.1.1.1 1995/10/24 10:21:52 cernlib
15 #include "geant321/pilot.h"
16 *CMZ : 3.21/04 23/02/95 14.53.54 by S.Giani
17 *-- Author : Christian Zeitnitz 21/07/92
19 C**************************************************************
24 C Purpose : setup cross-section tables and initialize pointer
29 C last modification: Changed in order to read new x-section file
32 C for details see MICAP manual ORNL/TM-10340
33 C*************************************************************
35 #include "geant321/mmicap.inc"
36 #include "geant321/mpoint.inc"
37 #include "geant321/minput.inc"
38 #include "geant321/mmass.inc"
39 #include "geant321/mconst.inc"
40 #include "geant321/cmagic.inc"
41 #include "geant321/cerrcm.inc"
42 #include "geant321/camass.inc"
44 #include "geant321/gccuts.inc"
45 #include "geant321/gcflag.inc"
46 #include "geant321/gcunit.inc"
47 C pointer to material/mixture bank (NMATE,JMATE)
48 #include "geant321/gcnum.inc"
50 COMMON / QUEST / IQUEST(100)
52 C Avogadro number multiplied by 1.E-24
53 PARAMETER(XNAVO = 0.60221367)
55 DIMENSION A(100),AGEA(100),Z(100),DEN(100),MID(100,2),IDI(20,2)
61 #if defined(CERNLIB_UNIX)||defined(CERNLIB_CRAY)||defined(CERNLIB_VAX)
64 LOGICAL OPENED,EXISTS,FMIST,FSINGL,FMIFL
65 C GEANT Particle IDs used to extract masses from GEANT
66 DATA IPID /14 , 13 , 8 , 7 , 9 , 5 , 6 , 45 , 46 , 49 , 47 , 1/
68 C Initialization flag of GEANT
69 C Loop for XMASS extracted from GCALOR routine CALINI K.L-P
73 CALL GFPART(IPID(I),NAP,ITR,AM,CH,TL,UB,NW)
77 C neutron energy cut (eV)
79 C get time cut off from GEANT
81 C temperature for thermal neutron xsection (Kelvin)
82 C only temporary constant !!!
88 C open MICAP I/O units
89 INQUIRE(UNIT=MICROS,OPENED=OPENED)
93 #if defined(CERNLIB_VAX)
95 INQUIRE(FILE=XSFILE,EXIST=EXISTS)
97 ISTAT = LIB$SYS_TRNLOG('CERN_ROOT',NALL,CHROOT,,,%VAL(0))
98 IF(ISTAT.EQ.1) XSFILE = 'CERN_ROOT:[LIB]xsneut95.dat'
100 INQUIRE(FILE=XSFILE,EXIST=EXISTS)
102 PRINT*,'**********************************'
103 PRINT*,'* G C A L O R *'
104 PRINT*,'* ----------- *'
105 PRINT*,'* File XSNEUT95.DAT not found *'
106 PRINT*,'* Program STOP *'
107 PRINT*,'* Check CERN_ROOT environment *'
108 PRINT*,'* variable *'
109 PRINT*,'**********************************'
112 OPEN(UNIT=MICROS,FILE=XSFILE, STATUS='OLD',READONLY)
114 #if defined(CERNLIB_UNIX)||defined(CERNLIB_CRAY)
115 XSFILE = 'xsneut95.dat'
116 INQUIRE(FILE=XSFILE,EXIST=EXISTS)
119 CALL GETENV('CERN_ROOT',CHROOT)
120 LNROOT = LNBLNK(CHROOT)
122 + XSFILE = CHROOT(1:LNROOT)//'/lib/xsneut95.dat'
124 INQUIRE(FILE=XSFILE,EXIST=EXISTS)
126 PRINT*,'**********************************'
127 PRINT*,'* G C A L O R *'
128 PRINT*,'* ----------- *'
129 PRINT*,'* File XSNEUT95.DAT not found *'
130 PRINT*,'* Program STOP *'
131 PRINT*,'* Check CERN_ROOT environment *'
132 PRINT*,'* variable *'
133 PRINT*,'**********************************'
136 OPEN(UNIT=MICROS,FILE=XSFILE,STATUS='OLD')
138 #if defined(CERNLIB_MSDOS)||defined(CERNLIB_WINNT)
140 CALL GETENVF('CERN_ROOT',CHROOT)
141 LNROOT = LNBLNK(CHROOT)
143 XSFILE = 'xsneut95.dat'
145 XSFILE = CHROOT(1:LNROOT)//'\\lib\\xsneut95.dat'
146 INQUIRE(FILE=XSFILE,EXIST=EXISTS)
147 IF(.NOT.EXISTS) XSFILE='xsneut95.dat'
149 OPEN(UNIT=MICROS,FILE=XSFILE)
151 #if defined(CERNLIB_IBMVM)
152 XSFILE = '\XSNEUT95 DAT *'
153 OPEN(UNIT=MICROS,FILE=XSFILE,STATUS='OLD')
156 C setup the link areas needed for x-section banks
157 CALL MZLINK(IXCONS,'MICTMP',LTEMP,LTEMP,LTEMP)
158 CALL MZLINK(IXCONS,'MMICAP',LMAG2,LMOX4,LMAG2)
159 CALL MZLINK(IXCONS,'MPOINT',LMAG1,LFP210,LMAG1)
164 C pointers into TEMP bank
167 NTMPNI = NTNAME + 1 + 80/4
169 NTDATS = NTCOMM + 1 + 80/4
170 NTLIST = NTDATS + 1 + 24/4
172 C read comment and date of xsection file
173 READ(NUNIT,'(A80,/,A24)') COMMEN,DATSTR
174 C read in material definition array
175 READ(NUNIT,'(I10)') NISO
176 NWW = NISO * 3 + 12 + NTLIST
177 C get temporary buffer
178 CALL CHKZEB(NWW,IXCONS)
180 C create a top level bank for the list of isotopes
181 CALL MZBOOK(IXCONS,LTEMP,LSUP,1,'TEMP',3,0,NWW,0,-1)
184 C create an additional bank in the linear structure TEMP
185 CALL MZBOOK(IXCONS,LT,LSUP,0,'TEMP',3,0,NWW,0,-1)
190 C store the unit number of the file in bank TEMP
191 IQ(LT + NTUNIT) = NUNIT
192 C store the file name in bank TEMP
193 CALL UCTOH(XSFILE,IQ(LT+NTNAME+1),4,LNBLNK(XSFILE))
194 IQ(LT + NTNAME) = LNBLNK(XSFILE)
195 C store the comment and date string in bank TEMP
196 IQ(LT + NTCOMM) = LNBLNK(COMMEN)
197 CALL UCTOH(COMMEN,IQ(LT+NTCOMM+1),4,LNBLNK(COMMEN))
198 IQ(LT + NTDATS) = LNBLNK(DATSTR)
199 CALL UCTOH(DATSTR,IQ(LT+NTDATS+1),4,LNBLNK(DATSTR))
201 LL = (I-1)*12 + LT + NTLIST
202 READ(NUNIT,'(12I6)') (IQ(L),L=LL,LL+11)
205 C get number of comment lines for different isotopes
206 READ(NUNIT,'(I10)') NCOM
209 CALL CHKZEB(NWW,IXCONS)
211 C create a top level bank for the isotope comments
212 CALL MZBOOK(IXCONS,LCISO,LCSUP,1,'CISO',3,0,NWW,0,-1)
215 C create an additional bank in the linear structure CISO
216 CALL MZBOOK(IXCONS,LC,LCSUP,0,'CISO',3,0,NWW,0,-1)
222 READ(NUNIT,'(I4,I4,A70)') IQ(LC+J),IQ(LC+J+1),
224 CALL UCTOH(CCOMM,IQ(LC+J+2),4,70)
227 C---------------------------------------------------------------------
228 C check the existence of secondary x-section files stored in bank MIFL
229 C real messy code !!! But its fortran after all !!! CZ Jan 95
231 IF(NUNIT.EQ.MICROS) THEN
234 IF(LMIFIL.GE.IQUEST(3) .AND. LMIFIL.LE.IQUEST(4)) THEN
235 CALL UHTOC(IQ(LMIFIL-4),4,CNAME,4)
236 IF(CNAME.EQ.'MIFL') FMIFL = .TRUE.
243 CALL UHTOC(IQ(IXSF+2),4,XSFILE,IQ(IXSF+1))
245 INQUIRE(FILE=XSFILE,EXIST=EXISTS)
248 PRINT*,' * MICAP : x-section file not found : ',XSFILE
252 C last name in the list ?
253 IF(IXSF-LMIFIL .GE. IQ(LMIFIL-1) ) FMIFL = .FALSE.
254 C find a free unit number (greater 31), and use it
256 INQUIRE(UNIT=I,OPENED=OPENED)
259 #if defined(CERNLIB_UNIX)||defined(CERNLIB_CRAY)||defined(CERNLIB_IBMVM)
260 OPEN(UNIT=I,FILE=XSFILE,STATUS='OLD')
262 #if defined(CERNLIB_VAX)
263 OPEN(UNIT=I,FILE=XSFILE,STATUS='OLD',READONLY)
269 PRINT *,'* MICAP : No more free units available !'
273 C---------------------------------------------------------------------
274 CALL VZERO(MATIDS,4000)
277 NUNIT = IQ(LT + NTUNIT)
282 IF(NIS.EQ.0) GOTO 100
283 IF(MATIDS(I,1,1).EQ.0) THEN
285 MATIDS(I,1,2) = NUNIT
286 C is the Z of the element correct?
287 ELSE IF(IQ(KK)/1000.EQ.I) THEN
288 C overwrite existing element with the one stored in new file
289 DO 70 J=2,MATIDS(I,1,1)+1
294 MATIDS(I,1,2) = NUNIT
300 C maximal 20 isotopes per element
303 MATIDS(I,J,1) = IQ(KK)
304 MATIDS(I,J,2) = IQ(KK+1)
312 C DEFINE CROSS SECTION DIMENSIONING VARIABLES
313 C NNR EQUALS THE NUMBER OF NEUTRON RECORDS
314 C NQ EQUALS THE NUMBER OF Q VALUES
315 C NGR EQUALS THE NUMBER OF GAMMA RECORDS
316 C SET THE DEFAULT VALUES FOR THE CURRENT CROSS SECTION DATA
321 C SET THE DEFAULT VALUES FOR THE NEUTRON, PROTON, DEUTERON,
322 C TRITON, HELIUM-3, AND ALPHA PARTICLE MASSES (IN EV)
329 C SET THE DEFAULT VALUES FOR THE NEUTRON, PROTON, DEUTERON,
330 C TRITON, HELIUM-3, AND ALPHA PARTICLE MASSES (IN AMU)
338 C now preprocess all materials xs
342 #if defined(CERNLIB_MDEBUG)
343 PRINT *,' MICAP-INI : setup materials '
344 PRINT '('' NMATE='',I20,'' JMATE='',I20)',NMATE,JMATE
345 PRINT '('' NTMED='',I20,'' JTMED='',I20)',NTMED,JTMED
347 C Check if material option bank MIST exists
350 IF(LMIST.GE.IQUEST(3) .AND. LMIST.LE.IQUEST(4)) THEN
351 CALL UHTOC(IQ(LMIST-4),4,CNAME,4)
352 IF(CNAME.EQ.'MIST') FMIST = .TRUE.
354 C 1. loop over tracking media -> get NMIX,MEDIA
357 IF(JTM.LE.0) GOTO 140
358 C valid tracking medium found get material number
359 C and get corresponding material parameters from JMATE structure
361 IF(IMA.LE.0 .OR. IMA.GT.NMATE) GOTO 140
362 C count number of elements and number of mixing operations
364 IF(JMA.LE.0) GOTO 140
365 IF(Q(JMA+6) .LE. 1.0 .OR. Q(JMA+6) .GE. 240.) GOTO 140
366 C Check if for material IMA single isotopes are selected
369 DO 110 KIM=1,IQ(LMIST-1),2
370 IF(IMA.EQ.IQ(LMIST+KIM).AND.IQ(LMIST+KIM+1).EQ.0) THEN
377 C get number of elements in material max = 100
378 KK = MIN(ABS(Q(JMA+11)),100.)
379 C relation between MICAP and GEANT material number
386 C check if more than one isotope has to taken into account for all
387 C elements in the mixture
389 IA = NINT(Q(JMIXT+K))
390 IZ = NINT(Q(JMIXT+K+KK))
391 CALL MATISO(IZ,IA,NNI,IDI,FSINGL,NUNIT)
397 CALL MATISO(IZ,IA,NNI,IDI,FSINGL,NUNIT)
402 C allocate ZEBRA bank for material information
403 NW = 9 * NMIX + MEDIA + 10
404 C define link area for MICAP banks in GCBANK
405 CALL CHKZEB(NW,IXCONS)
406 CALL MZBOOK(IXCONS,LMOMA,0,2,'MOME',0,0,NW,0,-1)
410 LFP10 = LGE2MO + MEDIA + 1
415 LFP140 = LFP14 + NMIX
416 LFP16 = LFP140 + NMIX
418 C 2. loop over tracking media
423 IF(JTM.LE.0) GOTO 230
424 C valid tracking medium found get material number
425 C and get corresponding material parameters from JMATE structure
427 #if defined(CERNLIB_MDEBUG)
428 PRINT '('' IMATE ='',I10)',IMA
430 IF(IMA.LE.0 .OR. IMA.GT.NMATE) GOTO 230
431 C count number of elements and number of mixing operations
433 IF(JMA.LE.0) GOTO 230
434 IF(Q(JMA+6) .LE. 1.0 .OR. Q(JMA+6) .GE. 240.) GOTO 230
435 C Check if for material IMA single isotopes are selected
438 DO 150 KIM=1,IQ(LMIST-1),2
439 IF(IMA.EQ.IQ(LMIST+KIM).AND.IQ(LMIST+KIM+1).EQ.0) THEN
447 C get number of elements in material max = 100
449 KK = MIN1(ABS(Q(JMA+11)),100.)
450 C relation between MICAP and GEANT material number
451 C check if medium IMA already stored (multiple tracking media)
454 IF(IQ(LGE2MO+KMI).EQ.IMA) THEN
458 CALL MATISO(IZ,IA,NNI,IDI,FSINGL,NUNIT)
463 IA = NINT(Q(JMIXT+K))
464 IZ = NINT(Q(JMIXT+K+KK))
465 CALL MATISO(IZ,IA,NNI,IDI,FSINGL,NUNIT)
474 IQ(LGE2MO+MEDIA1) = IMA
481 AMOL = Q(LQ(JMIXT-1) + 2)
482 XMOLCM = RHO1/AMOL*XNAVO
483 IA = NINT(Q(JMIXT+K))
484 IZ = NINT(Q(JMIXT+K+KK))
485 CALL MATISO(IZ,IA,NNI,IDI,FSINGL,NUNIT)
488 KKPOS = KPOS + KJ - 1
490 AGEA(KKPOS) = Q(JMIXT+K)
494 MID(KKPOS,1) = IDI(KJ,1)
497 IIA = IDI(KJ,1) - IIZ * 1000
498 IF(IIA.NE.0 .AND. NNI.GT.1.) THEN
499 A(KKPOS) = FLOAT(IIA)
501 A(KKPOS) = Q(JMIXT+K)
503 Z(KKPOS) = Q(JMIXT+K+KK)
504 WISO = FLOAT(IDI(KJ,2))/100.
505 WI = Q(JMIXT+K+2*KK)*AMOL/A(KKPOS)*WISO
506 DEN(KKPOS) = XMOLCM * WI
507 #if defined(CERNLIB_MDEBUG)
508 PRINT '('' MIXT: El #'',I3,'' A,Z :'',2F10.2,
509 + '' Rho='',F10.2,'' Den ='',F10.5)',
510 + KJ,A(KJ),Z(KJ),RHO1,DEN(KJ)
515 C element or compound
519 CALL MATISO(IZ,IA,NNI,IDI,FSINGL,NUNIT)
527 MID(KJ,1) = IDI(KJ,1)
530 IIA = IDI(KJ,1) - IIZ * 1000
531 IF(IIA.NE.0 .AND. NNI.GT.1.) THEN
537 WISO = FLOAT(IDI(KJ,2))/100.
538 DEN(KJ) = RHO1/A(KJ) * WISO *XNAVO
539 #if defined(CERNLIB_MDEBUG)
540 PRINT '('' ELEM: Iso #'',I3,'' A,Z :'',2F10.2,
541 + '' Rho='',F10.2,'' Den ='',F10.5)',
542 + KJ,A(KJ),Z(KJ),RHO1,DEN(KJ)
547 C fill MICAP material arrays
548 C actual number of isotopes given by KK2
550 DO 220 J = NMIX1 + 1, NMIX1 + KK2
551 IQ(LFP10+J-1) = MEDIA1
552 IQ(LFP11+J-1) = MID(J-NMIX1,1)
553 C check if bound hydrogen has been selected
554 IF(NINT(Z(J-NMIX1)).EQ.1.AND.KK.GT.1) IQ(LFP11+J-1) = 1000
555 Q(LFP12+J-1) = DEN(J-NMIX1)
556 IQ(LFP13+J-1) = NINT(Z(J-NMIX1))
557 Q(LFP14+J-1) = A(J-NMIX1)
558 Q(LFP140+J-1) = AGEA(J-NMIX1)
563 PRINT *,' GCALOR: NO tracking media found ===> STOP '
566 C read cross-sections and perform mixing and thinning
568 C close MICAP cross-section file(s)
571 CLOSE(UNIT=IQ(LT+NTUNIT))
574 C Drop temporary linear structures
575 CALL MZDROP(IXCONS,LTEMP,'L')
576 CALL MZDROP(IXCONS,LCISO,'L')