]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/gphys/gphxsi.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Some function moved to AliZDC
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gphxsi.F
CommitLineData
fe4da5cc 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
45c
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
304c 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
369c 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