Use fluka routine usrdci to get ion properties.
[u/mrichter/AliRoot.git] / TFluka / comscw_activity.f
CommitLineData
93d28366 1*$ CREATE COMSCW.FOR
2*COPY COMSCW
3*
4*=== comscw ===========================================================*
5*
6 DOUBLE PRECISION FUNCTION COMSCW ( IJ , XA , YA , ZA ,
7 & MREG , RULL , LLO , ICALL )
8
9 INCLUDE '(DBLPRC)'
10 INCLUDE '(DIMPAR)'
11 INCLUDE '(IOUNIT)'
12* *
13*----------------------------------------------------------------------*
14* *
15* Special multiplication factors for activity scoring. *
16* *
17* Note: Data files 'LE-CH.dat' and 'LE-EC.dat' read on unit 20. *
18* *
19* *
20* SDUM = 'abcdefgh' Lesco *
21* *
22* 'ab' = 'LE' division by exemption limits activated *
23* *
24* 'c' = 'C' : division by CERN zonage limits *
25* 'def' = '10 ' : division by 1/10 of the value 10 *
26* 'def' = '100' : division by 1/100 of the value 100 *
27* otherwise: division by the value itself 1 *
28* 'c' = 'S' : division by Swiss LE values 2 *
29* 'c' = 'E' : division by European concentration limits 3 *
30* *
31* Comscw = 1. for selected isotope *
32* = 0. otherwise *
33* *
34* 'ab' = 'IS' selection of individual isotops activated 9999 *
35* *
36* 'cd' = symbol of isotope *
37* 'efg' = mass number of isotope *
38* 'h' = 'm' : metastable state *
39* *
40* Comscw = 1. for selected isotope *
41* = 0. otherwise *
42* *
43* otherwise Comscw = 1. -1 *
44* *
45* Example(s): *
46* *
47* SDUM = 'LEC10xxx' : division by 1/10 of LE(EU) or by Swiss LE *
48* = 'LESxxxxx' : division by Swiss LE *
49* = 'ISV 48x' : scoring of activity for V48 *
50* = 'ISMn 52m' : scoring of activity for mMn52 *
51* *
52* Note: 'x' refers to characters without specific meaning *
53* *
54* Input variables: *
55* *
56* Ij = (generalized) particle code *
57* Xa,Ya,Za = position *
58* Mreg = region number *
59* Rull = amount to be deposited *
60* Llo = particle generation *
61* Icall = call id *
62* *
63* Output variables: *
64* *
65* Comscw = factor the scored amount will be multiplied by *
66* Lsczer = logical flag, if true no amount will be scored *
67* regardless of Comscw *
68* *
69* Useful variables (common SCOHLP): *
70* *
71* Energy/Star binnings/scorings (Comscw): *
72* ISCRNG = 1 --> Energy density binning *
73* ISCRNG = 2 --> Star density binning *
74* ISCRNG = 3 --> Residual nuclei scoring *
75* ISCRNG = 4 --> Momentum transfer density binning *
76* ISCRNG = 5 --> Activity density binning *
77* JSCRNG = # of the binning *
78* *
79* Useful variables (common SOUEVT): *
80* *
81* X,Y,Zsoevt(i) = position of the i_th source particle *
82* TX,Y,Zsoev(i) = direction of the i_th source particle *
83* Wtsoev(i) = weight of the i_th source particle *
84* Pmsoev(i) = momentum of the i_th source particle *
85* Tksoev(i) = kin. energy of the i_th source particle *
86* Agsoev(i) = age of the i_th source particle *
87* Aksoev(i) = Kaon ampl. of the i_th source particle *
88* Ussoev(i) = user var. of the i_th source particle *
89* Ijsoev(i) = identity of the i_th source particle *
90* Nrsoev(i) = region of the i_th source particle *
91* Nlsoev(i) = lattice of the i_th source particle *
92* Npsoev = number of the source particles *
93*----------------------------------------------------------------------*
94*
95 INCLUDE '(FLKMAT)'
96 INCLUDE '(SCOHLP)'
97 INCLUDE '(SOUEVT)'
98*
99 INCLUDE '(RSNCCM)'
100 INCLUDE '(USRBIN)'
101 DIMENSION IZSCO(MXUSBN),IASCO(MXUSBN),ISSCO(MXUSBN),
102 & LESCO(MXUSBN)
103*
104 CHARACTER CISOIN*2,CISO*2,CSET*10,CDUM*4,CA*3
105 PARAMETER (IZMAX = 109,
106 & IAZMIN = 2,
107 & IAZMAX = 160)
108 DIMENSION CISO(IZMAX)
109 DIMENSION XLESWS(IZMAX,IAZMIN:IAZMAX,2),
110 & XLEECO(IZMAX,IAZMIN:IAZMAX,2),
111 & XLEO10(IZMAX,IAZMIN:IAZMAX,2)
112 DATA CISO /
113 & 'H ','He','Li','Be','B ','C ','N ','O ','F ','Ne','Na',
114 & 'Mg','Al','Si','P ','S ','Cl','Ar','K ','Ca','Sc','Ti',
115 & 'V ','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As',
116 & 'Se','Br','Kr','Rb','Sr','Y ','Zr','Nb','Mo','Tc','Ru',
117 & 'Rh','Pd','Ag','Cd','In','Sn','Sb','Te','I ','Xe','Cs',
118 & 'Ba','La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy',
119 & 'Ho','Er','Tm','Yb','Lu','Hf','Ta','W ','Re','Os','Ir',
120 & 'Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn','Fr','Ra',
121 & 'Ac','Th','Pa','U ','Np','Pu','Am','Cm','Bk','Cf','Es',
122 & 'Fm','Md','No','Lr','Rf','Ha','Sg','Ns','Hs','Mt'/
123*
124 LOGICAL LFIRST,LFOUND
125 DATA LFIRST /.TRUE./
126*
127 LSCZER = .FALSE.
128 COMSCW = ONEONE
129
130 IF (LFIRST) THEN
131
132 WRITE(LUNOUT,'(A)') ' COMSCW: activity weighting activated'
133 LFIRST = .FALSE.
134*
135*----------------------------------------------------------------------*
136* Initialization
137*
138* LE data taken from Appendix 3 column 9,
139* Ordonnance sur la radioprotection (ORaP) du 22 juin 1994
140* (etat au 4 avril 2000)
141*
142* (data file in Bq/kg!)
143*
144 DO 104 IZ=1,IZMAX
145 DO 105 IAZ=IAZMIN,IAZMAX
146 XLESWS(IZ,IAZ,1) = ZERZER
147 XLESWS(IZ,IAZ,2) = ZERZER
148 105 CONTINUE
149 104 CONTINUE
150*
151 OPEN(20,FILE='LE-CH.dat',STATUS='UNKNOWN')
152 NISO = 0
153 IZ0 = 1
154 100 CONTINUE
155 READ(20,*,END=101) CISOIN,IA,XLIMIT
156* convert from Bq/kg into Bq/g
157 XLIMIT = 1.0D-3*XLIMIT
158 LFOUND = .FALSE.
159 DO 102 IZ=IZ0,IZMAX
160 IF (CISOIN.EQ.CISO(IZ)) THEN
161 IF (IA.LT.0) THEN
162 IS = 2
163 IA = ABS(IA)
164 ELSE
165 IS = 1
166 ENDIF
167 IAZ = IA-IZ
168 IF ((IAZ.LT.IAZMIN).OR.(IAZ.GT.IAZMAX)) THEN
169 WRITE(LUNOUT,*)
170 & ' COMSCW: warning! Iaz out of allowed range: ',
171 & IAZ, IAZMIN,IAZMAX
172 STOP
173 ENDIF
174 IF (XLESWS(IZ,IAZ,IS).GT.ZERZER) THEN
175 WRITE(LUNOUT,*)
176 & ' COMSCW: warning! two entries for this isotope: ',
177 & CISOIN,IA,XLIMIT
178 STOP
179 ELSE
180 XLESWS(IZ,IAZ,IS) = XLIMIT
181 ENDIF
182 NISO = NISO+1
183 IZ0 = IZ
184 LFOUND = .TRUE.
185 GOTO 103
186 ENDIF
187 102 CONTINUE
188 WRITE(LUNOUT,*)
189 & ' COMSCW: isotope not recognized: ',CISOIN,IA,XLIMIT
190 STOP
191 103 CONTINUE
192 GOTO 100
193 101 CONTINUE
194 CLOSE(20)
195*
196* LE data taken from
197* Official journal of the European Communities L159, 29 June 1996
198* Council Directive 96/29/Euratom
199*
200* (data file in Bq/g!)
201*
202 DO 204 IZ=1,IZMAX
203 DO 205 IAZ=IAZMIN,IAZMAX
204 XLEECO(IZ,IAZ,1) = ZERZER
205 XLEECO(IZ,IAZ,2) = ZERZER
206 205 CONTINUE
207 204 CONTINUE
208*
209 OPEN(20,FILE='LE-EC.dat',STATUS='UNKNOWN')
210 NISO = 0
211 IZ0 = 1
212 200 CONTINUE
213 READ(20,*,END=201) CISOIN,IA,XLIMIT,IFLAG
214 LFOUND = .FALSE.
215 DO 202 IZ=IZ0,IZMAX
216 IF (CISOIN.EQ.CISO(IZ)) THEN
217 IF (IA.LT.0) THEN
218 IS = 2
219 IA = ABS(IA)
220 ELSE
221 IS = 1
222 ENDIF
223 IAZ = IA-IZ
224 IF ((IAZ.LT.IAZMIN).OR.(IAZ.GT.IAZMAX)) THEN
225 WRITE(LUNOUT,*)
226 & ' COMSCW: warning! Iaz out of allowed range: ',
227 & IAZ, IAZMIN,IAZMAX
228 STOP
229 ENDIF
230 IF (XLEECO(IZ,IAZ,IS).GT.ZERZER) THEN
231 WRITE(LUNOUT,*)
232 & ' COMSCW: warning! two entries for this isotope: ',
233 & CISOIN,IA,XLIMIT
234 STOP
235 ELSE
236 XLEECO(IZ,IAZ,IS) = XLIMIT
237* zero entries with Swiss values
238 IF (IFLAG.EQ.1) XLEECO(IZ,IAZ,IS) = ZERZER
239 ENDIF
240 NISO = NISO+1
241 IZ0 = IZ
242 LFOUND = .TRUE.
243 GOTO 203
244 ENDIF
245 202 CONTINUE
246 WRITE(LUNOUT,*)
247 & ' COMSCW: isotope not recognized: ',CISOIN,IA,XLIMIT
248 STOP
249 203 CONTINUE
250 GOTO 200
251 201 CONTINUE
252 CLOSE(20)
253 DO 404 IZ=1,IZMAX
254 DO 405 IAZ=IAZMIN,IAZMAX
255 DO 406 IS=1,2
256 XLEO10(IZ,IAZ,IS) = ZERZER
257 IF (XLEECO(IZ,IAZ,IS).GT.ZERZER) THEN
258 XLEO10(IZ,IAZ,IS) = XLEECO(IZ,IAZ,IS)/10.0D0
259 ELSE
260 XLEO10(IZ,IAZ,IS) = XLESWS(IZ,IAZ,IS)
261 ENDIF
262 406 CONTINUE
263 405 CONTINUE
264 404 CONTINUE
265*
266 DO 500 I=1,MXUSBN
267 IZSCO(I) = 0
268 IASCO(I) = 0
269 ISSCO(I) = 0
270 LESCO(I) = 0
271 500 CONTINUE
272
273 ENDIF
274*
275*----------------------------------------------------------------------*
276* Online weighting
277*
278 IF ( ISCRNG .EQ. 5 ) THEN
279*
280* determine type of weighting
281 IF (LESCO(JSCRNG).EQ.0) THEN
282 CSET = TITUSB(JSCRNG)
283 IF (CSET(1:2).EQ.'IS') THEN
284 LESCO(JSCRNG) = 9999
285 DO 700 IZ=1,IZMAX
286 IF (CSET(3:4).EQ.CISO(IZ)) THEN
287 IZSCO(JSCRNG) = IZ
288 READ(CSET,'(A4,A3)') CDUM,CA
289 READ(CA,'(I3)') IASCO(JSCRNG)
290 IF (CSET(8:8).EQ.'m') THEN
291 ISSCO(JSCRNG) = 2
292 ELSE
293 ISSCO(JSCRNG) = 1
294 ENDIF
295 ENDIF
296 700 CONTINUE
297 IF ((IZSCO(JSCRNG).LE.0).OR.(IASCO(JSCRNG).LE.0)
298 & .OR.(ISSCO(JSCRNG).LE.0)) THEN
299 WRITE(LUNOUT,*) ' COMSCW: unknown isotope, Z,A,S = ',
300 & IZSCO(JSCRNG),IASCO(JSCRNG),ISSCO(JSCRNG)
301 STOP
302 ENDIF
303 ELSEIF(CSET(1:2).EQ.'LE') THEN
304 IF (CSET(3:3).EQ.'C') THEN
305 IF (CSET(4:6).EQ.'100') THEN
306 LESCO(JSCRNG) = 100
307 ELSEIF (CSET(4:6).EQ.'10') THEN
308 LESCO(JSCRNG) = 10
309 ELSE
310 LESCO(JSCRNG) = 1
311 ENDIF
312 ELSEIF (CSET(3:3).EQ.'S') THEN
313 LESCO(JSCRNG) = 2
314 ELSEIF (CSET(3:3).EQ.'E') THEN
315 LESCO(JSCRNG) = 3
316 ELSE
317 WRITE(LUNOUT,*) ' COMSCW: unknown LE set ',CSET(3:3)
318 STOP
319 ENDIF
320 ELSE
321 LESCO(JSCRNG) = -1
322 ENDIF
323 WRITE(LUNOUT,1000) ' COMSCW: scoring ',JSCRNG,CSET,
324 & ' weighted with properties ',
325 & LESCO(JSCRNG),IZSCO(JSCRNG),IASCO(JSCRNG),ISSCO(JSCRNG)
326 1000 FORMAT(A,I3,2A,I5,3I4)
327 ENDIF
328*
329* obtain present isotope from common block
330 JA = IARSDL(1)
331 JZ = IZRSDL(1)
332 JS = ISRSDL(1)+1
333 IF (JS.GT.2) JS = 2
334 JAZ = JA-JZ
335*
336 COMSCW = ZERZER
337 IF (LESCO(JSCRNG).EQ.-1) THEN
338 COMSCW = ONEONE
339 ELSEIF (LESCO(JSCRNG).EQ.9999) THEN
340 IF ((JA.EQ.IASCO(JSCRNG)).AND.(JZ.EQ.IZSCO(JSCRNG)).AND.
341 & (JS.EQ.ISSCO(JSCRNG))) THEN
342 COMSCW = ONEONE
343 ELSE
344 COMSCW = ZERZER
345 ENDIF
346 ELSEIF ((LESCO(JSCRNG).EQ.1).OR.(LESCO(JSCRNG).EQ.10).OR.
347 & (LESCO(JSCRNG).EQ.100)) THEN
348 FACT = 10.0D0/DBLE(LESCO(JSCRNG))
349 IF (XLEO10(JZ,JAZ,JS).GT.ZERZER)
350 & COMSCW = ONEONE/(FACT*XLEO10(JZ,JAZ,JS))
351 ELSEIF (LESCO(JSCRNG).EQ.2) THEN
352 IF (XLESWS(JZ,JAZ,JS).GT.ZERZER)
353 & COMSCW = ONEONE/XLESWS(JZ,JAZ,JS)
354 ELSEIF (LESCO(JSCRNG).EQ.3) THEN
355 IF (XLEECO(JZ,JAZ,JS).GT.ZERZER)
356 & COMSCW = ONEONE/XLEECO(JZ,JAZ,JS)
357 ELSE
358 WRITE(LUNOUT,*) ' COMSCW: invalid option ',LESCO(JSCRNG)
359 STOP
360 ENDIF
361
362 ENDIF
363
364 RETURN
365*=== End of function Comscw ===========================================*
366 END