This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / neutron / photon.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:58  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/04 23/02/95  14.46.01  by  S.Giani
11 *-- Author :
12       SUBROUTINE PHOTON(D,LD,IDICTS,LDICT,NTX,NTS,IGCBS,LGCB,
13      +      AWR,IGCBS2,LGCB2,LR,IGAMS,LGAM,QI,ID,IIN,LRI,SIGN)
14 C       THIS ROUTINE CONTROLS THE GENERATION AND STORAGE OF ALL
15 C       PHOTONS PRODUCED BY THE NEUTRON INTERACTIONS.  WHERE DATA
16 C       PERMITS, THE PHOTON PRODUCED IS DIRECTLY COUPLED TO THE
17 C       NEUTRON REACTION OCCURING.
18 #include "geant321/minput.inc"
19 #include "geant321/mconst.inc"
20 #include "geant321/mnutrn.inc"
21 #include "geant321/mapoll.inc"
22 #include "geant321/mcross.inc"
23 #include "geant321/mpstor.inc"
24 #include "geant321/mmicab.inc"
25       DIMENSION IDICTS(NNR,NNUC),LDICT(NNR,NNUC),NTX(*),NTS(*),
26      +  IGCBS(NGR,NNUC),LGCB(NGR,NNUC),AWR(*),IGCBS2(NGR,NNUC),
27      +  LGCB2(NGR,NNUC),LR(NQ,NNUC),IGAMS(*),LGAM(*),D(*),LD(*)
28       SAVE
29 C flag to mark call to SECEGY = 1 or PARTXS = 2 for EP   CZ 13/8/92
30       IEP = 0
31 C       INITIALIZE THE PHOTON ENERGY TO ZERO IN CASE NO PHOTON IS
32 C       CHOSEN (THIS IS NECESSARY BECAUSE OF ENDF INCONSISTENCY)
33       EG=0.0
34 C       INITIALIZE THE PARAMETERS USED IN THE SELECTION PROCESS
35       MT=0
36       IMT=0
37       NUMBG=0
38       XSIG2=0.0
39       XSIG=0.0
40       SIGMT3=0.0
41       SIGP=0.0
42       AWRI=AWR(IIN)
43       NNTX=NTX(IIN)
44       NNTS=NTS(IIN)
45       L=2*NNTX+2*NNTS
46 C       NO PHOTON DATA PRESENT (IF L=0)
47       IF(L.EQ.0)GO TO 360
48       LX=2*NNTX
49       LS=LX+1
50 C       DETERMINE THE NEUTRON REACTION MT NUMBER
51       IF(ID.EQ.8)MT=16
52       IF(ID.EQ.9)MT=17
53       IF(ID.EQ.10)MT=18
54       IF(ID.EQ.11)MT=22
55       IF(ID.EQ.12)MT=24
56       IF(ID.EQ.13)MT=28
57       IF((ID.GE.14).AND.(ID.LE.54))MT=51
58       IF(ID.EQ.55)MT=102
59       IF(ID.EQ.56)MT=103
60       IF(ID.EQ.57)MT=104
61       IF(ID.EQ.58)MT=105
62       IF(ID.EQ.59)MT=106
63       IF(ID.EQ.60)MT=107
64       IF(ID.EQ.61)MT=108
65       IF(ID.EQ.62)MT=109
66       IF(ID.EQ.63)MT=111
67       IF(ID.EQ.64)MT=112
68       IF(ID.EQ.65)MT=113
69       IF(ID.EQ.66)MT=114
70 C       DETERMINE WHICH DISCRETE INELASTIC SCATTERING LEVEL OCCURRED
71       IF(MT.NE.51)GO TO 130
72       IMT=ID-14
73       MT=MT+IMT
74 C       RESET THE MT NUMBER IF AN LR-FLAG IS INVOLKED
75       IF(LRI.EQ.22)MT=22
76       IF(LRI.EQ.23)MT=23
77       IF(LRI.EQ.28)MT=28
78 C       CHECK PHOTON PRODUCTION DICTIONARY TO SEE IF THERE IS PHOTON
79 C       DATA CORRESPONDING TO THE NEUTRON MT REACTION THAT OCCURRED
80       DO 10 IX=1,NNTX
81          MTG=LGCB(2*IX-1,IIN)
82          IF(MTG.EQ.MT)GO TO 30
83    10 CONTINUE
84    20 IF(LRI.EQ.22)GO TO 190
85       IF(LRI.EQ.23)GO TO 190
86       IF(LRI.EQ.28)GO TO 190
87       GO TO 70
88 C       PHOTON DATA FOUND CORRESPONDING TO NEUTRON MT REACTION
89    30 L1=LGCB2(2*IX,IIN)
90       IF(L1.EQ.0)GO TO 370
91       LS1=IGCBS2(2*IX,IIN)+LMOX4
92       LEN=L1/2
93       CALL TBSPLT(D(LS1),EOLD,LEN,SIGP)
94       IF(SIGP.EQ.0.0)GO TO 190
95       LS2=IGCBS(2*IX,IIN)+LMOX2
96 C       DETERMINE EXIT PHOTON ENERGY (EP)
97       CALL PARTXS(D(LS2),D(LS2),EOLD,SIGP,EP)
98       IEP = 2
99       IF(EP.GT.0.0)GO TO 60
100 C       DISCRETE PHOTON ENERGY WAS NOT SELECTED (EP=0.0)
101 C       CHECK SECONDARY PHOTON DISTRIBUTION (FILE 15) FOR EP
102       DO 40 IS=1,NNTS
103          MTGS=LGCB(LX+2*IS-1,IIN)
104          IF(MTGS.EQ.MT)GO TO 50
105    40 CONTINUE
106 C no file 15 found and EP=0 in PARTXS -> try MT=4 etc
107       GO TO 20
108    50 L1=LGCB(LX+2*IS,IIN)
109       IF(L1.EQ.0)GO TO 380
110       LS3=IGCBS(LX+2*IS,IIN)+LMOX2
111 C       DETERMINE EXIT PHOTON ENERGY (EP)
112       CALL SECEGY(EP,D(LS3),EOLD,D(LS3))
113       IEP = 1
114 C       DETERMINE THE PHOTON MULTIPLICITY (YP)
115 C       RECALCULATE THE DENOMINATOR USED IN CALCULATING THE
116 C       PHOTON MULTIPLICITY TO ACCOUNT FOR THE LR-FLAGS
117    60 IF(LRI.EQ.22)CALL LRNORM(D,D,IDICTS,LDICT,LR,EOLD,MT,IIN,SIGN)
118       IF(LRI.EQ.23)CALL LRNORM(D,D,IDICTS,LDICT,LR,EOLD,MT,IIN,SIGN)
119       IF(LRI.EQ.28)CALL LRNORM(D,D,IDICTS,LDICT,LR,EOLD,MT,IIN,SIGN)
120       YP=SIGP/SIGN
121       GO TO 330
122 C       THE DISCRETE INELASTIC LEVEL PHOTON DATA WAS NOT FOUND
123 C       CHECK THE PHOTON PRODUCTION DICTIONARY TO SEE IF THERE IS
124 C       PHOTON DATA CORRESPONDING TO MT=4
125    70 DO 80 IX=1,NNTX
126          MTG=LGCB(2*IX-1,IIN)
127          IF(MTG.EQ.4)GO TO 90
128    80 CONTINUE
129       GO TO 190
130 C       PHOTON DATA FOUND CORRESPONDING TO MT=4
131    90 L1=LGCB2(2*IX,IIN)
132       IF(L1.EQ.0)GO TO 370
133       LS1=IGCBS2(2*IX,IIN)+LMOX4
134       LEN=L1/2
135       CALL TBSPLT(D(LS1),EOLD,LEN,SIGP)
136       IF(SIGP.EQ.0.0)GO TO 190
137       LS2=IGCBS(2*IX,IIN)+LMOX2
138 C       DETERMINE EXIT PHOTON ENERGY (EP)
139       CALL PARTXS(D(LS2),D(LS2),EOLD,SIGP,EP)
140       IEP = 2
141       IF(EP.GT.0.0)GO TO 120
142 C       DISCRETE PHOTON ENERGY WAS NOT SELECTED (EP=0.0)
143 C       CHECK SECONDARY PHOTON DISTRIBUTION (FILE 15) FOR EP
144       DO 100 IS=1,NNTS
145          MTGS=LGCB(LX+2*IS-1,IIN)
146          IF(MTGS.EQ.4)GO TO 110
147   100 CONTINUE
148       GO TO 380
149   110 L1=LGCB(LX+2*IS,IIN)
150       IF(L1.EQ.0)GO TO 380
151       LS3=IGCBS(LX+2*IS,IIN)+LMOX2
152 C       DETERMINE EXIT PHOTON ENERGY (EP)
153       CALL SECEGY(EP,D(LS3),EOLD,D(LS3))
154       IEP = 1
155 C       DETERMINE THE PHOTON MULTIPLICITY (YP)
156 C       RECALCULATE THE DENOMINATOR USED IN CALCULATING THE
157 C       PHOTON MULTIPLICITY TO ACCOUNT FOR THE LR-FLAGS
158   120 MT=4
159       CALL LRNORM(D,D,IDICTS,LDICT,LR,EOLD,MT,IIN,SIGNIS)
160       YP=SIGP/SIGNIS
161       GO TO 330
162 C       CHECK PHOTON PRODUCTION DICTIONARY TO SEE IF THERE IS PHOTON
163 C       DATA CORRESPONDING TO THE NEUTRON MT REACTION THAT OCCURRED
164   130 DO 140 IX=1,NNTX
165          MTG=LGCB(2*IX-1,IIN)
166          IF(MTG.EQ.MT)GO TO 150
167   140 CONTINUE
168       GO TO 190
169 C       PHOTON DATA FOUND CORRESPONDING TO NEUTRON MT REACTION
170   150 L1=LGCB2(2*IX,IIN)
171       IF(L1.EQ.0)GO TO 370
172       LS1=IGCBS2(2*IX,IIN)+LMOX4
173       LEN=L1/2
174       CALL TBSPLT(D(LS1),EOLD,LEN,SIGP)
175       IF(SIGP.EQ.0.0)GO TO 190
176       LS2=IGCBS(2*IX,IIN)+LMOX2
177 C       DETERMINE EXIT PHOTON ENERGY (EP)
178       CALL PARTXS(D(LS2),D(LS2),EOLD,SIGP,EP)
179       IEP = 2
180       IF(EP.GT.0.0)GO TO 180
181 C       DISCRETE PHOTON ENERGY WAS NOT SELECTED (EP=0.0)
182 C       CHECK SECONDARY PHOTON DISTRIBUTION (FILE 15) FOR EP
183       DO 160 IS=1,NNTS
184          MTGS=LGCB(LX+2*IS-1,IIN)
185          IF(MTGS.EQ.MT)GO TO 170
186   160 CONTINUE
187       GO TO 380
188   170 L1=LGCB(LX+2*IS,IIN)
189       IF(L1.EQ.0)GO TO 380
190       LS3=IGCBS(LX+2*IS,IIN)+LMOX2
191 C       DETERMINE EXIT PHOTON ENERGY (EP)
192       CALL SECEGY(EP,D(LS3),EOLD,D(LS3))
193       IEP = 1
194 C       DETERMINE THE PHOTON MULTIPLICITY (YP)
195   180 YP=SIGP/SIGN
196       GO TO 330
197 C       NO PHOTON DATA WAS FOUND FOR THE PARTICULAR NEUTRON MT
198 C       REACTION OR FOR NEUTRON MT=4, THEREFORE CHECK THE PHOTON
199 C       PRODUCTION DICTIONARY TO SEE IF THERE IS PHOTON DATA
200 C       CORRESPONDING TO MT=3 (THE CATCH-ALL MT)
201   190 DO 200 IX=1,NNTX
202          MTG=LGCB(2*IX-1,IIN)
203          IF(MTG.EQ.3)GO TO 210
204   200 CONTINUE
205       GO TO 360
206 C       PHOTON DATA FOUND CORRESPONDING TO MT=3
207   210 L1=LGCB2(2*IX,IIN)
208       IF(L1.EQ.0)GO TO 370
209       LS1=IGCBS2(2*IX,IIN)+LMOX4
210       LEN=L1/2
211       CALL TBSPLT(D(LS1),EOLD,LEN,SIGP)
212       IF(SIGP.EQ.0.0)GO TO 360
213       LS2=IGCBS(2*IX,IIN)+LMOX2
214 C       DETERMINE EXIT PHOTON ENERGY (EP)
215       CALL PARTXS(D(LS2),D(LS2),EOLD,SIGP,EP)
216       IEP = 2
217       IF(EP.GT.0.0)GO TO 240
218 C       DISCRETE PHOTON ENERGY WAS NOT SELECTED (EP=0.0)
219 C       CHECK SECONDARY PHOTON DISTRIBUTION (FILE 15) FOR EP
220       DO 220 IS=1,NNTS
221          MTGS=LGCB(LX+2*IS-1,IIN)
222          IF(MTGS.EQ.3)GO TO 230
223   220 CONTINUE
224       GO TO 380
225   230 L1=LGCB(LX+2*IS,IIN)
226       IF(L1.EQ.0)GO TO 380
227       LS3=IGCBS(LX+2*IS,IIN)+LMOX2
228 C       DETERMINE EXIT PHOTON ENERGY (EP)
229       CALL SECEGY(EP,D(LS3),EOLD,D(LS3))
230       IEP = 1
231 C       THE PHOTON WAS SELECTED FROM PHOTON DATA FOR MT=3
232 C       TO OBTAIN THE CORRECT MULTIPLICITY, THE NEUTRON CROSS
233 C       SECTION FOR MT=3 MUST BE ADJUSTED TO REPRESENT THE SAME
234 C       DATA AS MT=3 DOES IN THE PHOTON DATA
235   240 ID=2
236 C       OBTAIN NEUTRON ELASTIC SCATTERING CROSS SECTION
237       L1=LDICT(ID,IIN)
238       IF(L1.EQ.0)GO TO 250
239       LS1=IDICTS(ID,IIN)+LMOX2
240       LEN=L1/2
241       CALL TBSPLT(D(LS1),EOLD,LEN,XSIG2)
242 C       SUBTRACT THE ELASTIC SCATTERING CROSS SECTION FROM THE TOTAL
243 C       CROSS SECTION TO OBTAIN BASE NEUTRON MT=3 REACTION
244       SIGMT3=SIGT-XSIG2
245       GO TO 260
246   250 SIGMT3=SIGT
247   260 CONTINUE
248 C       SCAN THE PHOTON PRODUCTION DICTIONARY FOR ALL MT NUMBERS
249 C       NOT EQUAL TO MT=3
250       DO 300 IX=1,NNTX
251          MTG=LGCB(2*IX-1,IIN)
252          IF(MTG.EQ.3)GO TO 300
253          L1=LGCB2(2*IX,IIN)
254          IF(L1.EQ.0)GO TO 370
255          LS1=IGCBS2(2*IX,IIN)+LMOX4
256          LEN=L1/2
257          CALL TBSPLT(D(LS1),EOLD,LEN,SIGEX)
258 C     IF THE TOTAL PHOTON PRODUCTION CROSS SECTION IS ZERO AT
259 C     THE NEUTRON ENERGY, THEN THE NEUTRON CROSS SECTION SHOULD
260 C     NOT BE SUBTRACTED FROM MT3 TO MAINTAIN PROPER NORMALIZATION
261          IF(SIGEX.EQ.0.0)GO TO 300
262 C     SET THE NEUTRON DICTIONARY ID NUMBER CORRESPONDING TO MTG
263          IF((MTG.LT.51).OR.(MTG.GT.91))GO TO 270
264          ID=14
265          IMT3=MTG-51
266          ID=ID+IMT3
267   270    IF(MTG.EQ.4)ID=3
268          IF(MTG.EQ.16)ID=8
269          IF(MTG.EQ.17)ID=9
270          IF(MTG.EQ.18)ID=10
271          IF(MTG.EQ.22)ID=11
272          IF(MTG.EQ.24)ID=12
273          IF(MTG.EQ.28)ID=13
274          IF(MTG.EQ.102)ID=55
275          IF(MTG.EQ.103)ID=56
276          IF(MTG.EQ.104)ID=57
277          IF(MTG.EQ.105)ID=58
278          IF(MTG.EQ.106)ID=59
279          IF(MTG.EQ.107)ID=60
280          IF(MTG.EQ.108)ID=61
281          IF(MTG.EQ.109)ID=62
282          IF(MTG.EQ.111)ID=63
283          IF(MTG.EQ.112)ID=64
284          IF(MTG.EQ.113)ID=65
285          IF(MTG.EQ.114)ID=66
286 C     OBTAIN THE NEUTRON CROSS SECTION CORRESPONDING TO MTG AND
287 C     SUBTRACT IT OFF OF THE BASE NEUTRON MT=3 CROSS SECTION
288          L1=LDICT(ID,IIN)
289          IF(L1.EQ.0)GO TO 280
290          LS1=IDICTS(ID,IIN)+LMOX2
291          LEN=L1/2
292          CALL TBSPLT(D(LS1),EOLD,LEN,XSIG)
293          GO TO 290
294   280    XSIG=0.0
295   290    SIGMT3=SIGMT3-XSIG
296          IF(SIGMT3.LE.0.0)GO TO 310
297   300 CONTINUE
298 C    DETERMINE THE PHOTON MULTIPLICITY (YP)
299       YP=SIGP/SIGMT3
300       IF(YP.GE.100.0)GO TO 310
301       GO TO 330
302   310 CONTINUE
303 C       THIS SECTION OF CODING IS INCLUDED TO ACCOUNT FOR ANY
304 C       ENDF/B DATA INCONSISTENCY WHICH COULD YIELD A PHOTON OF
305 C       CONSIDERABLE WEIGHT.  THE FOLLOWING CODING WILL SAMPLE THE
306 C       PHOTON WEIGHT FROM THE GENERAL PHOTON YIELD ARRAY AND
307 C       ADJUST THE WEIGHT TO PHOTONS PER NON-ELASTIC COLLISION.
308       L1=LGAM(IIN)
309       IF(L1.EQ.0)GO TO 320
310       LS1=IGAMS(IIN)+LMOX2
311       LEN=L1/2
312       CALL TBSPLT(D(LS1),EOLD,LEN,YP)
313       YP=(YP*SIGT)/(SIGT-XSIG2)
314       GO TO 330
315   320 YP=1.00
316 C       THE FOLLOWING SECTION OF CODING IS INCLUDED TO DISTRIBUTE
317 C       THE WEIGHT ENDF/B-V DATA MAY GIVE A PARTICULAR PHOTON.
318 C       FOR EXAMPLE, ENDF/B-V DATA MAY ASSIGN A MULITPLICITY OF
319 C       75 TO A PARTICULAR PHOTON.  BECAUSE SUCH A PHOTON COULD
320 C       CONSIDERABLY MODIFY THE RESULTS OF A DETECTOR RESPONSE, THE
321 C       MULTIPLICITY (PHOTON WEIGHT) IS DISTRIBUTED TO SEVERAL
322 C       PHOTONS (SPLITTING OF SORTS) WITH BOTH WEIGHT AND ENERGY
323 C       BEING CONSERVED.  THIS RARELY OCCURS BUT IS NECESSARY.
324   330 CONTINUE
325 C poisson distributed photon multiplicity CZ 13.8.92
326       IGTRY=0
327       MGPAR=INT(FLOAT(MAXPAR)*0.7)
328   340 CALL GPOISS(YP,NUMBG,1)
329       IGTRY=IGTRY+1
330       IF(NUMBG.GT.INT(4.*YP).OR.
331      +   NUMBG.GT.MGPAR.AND.IGTRY.LT.5) GOTO 340
332       NUMBG=MIN(NUMBG,MGPAR)
333 C Allow 0 Photond to be generated
334       IF(NUMBG.EQ.0) RETURN
335       EPTOT = YP*EP
336       EPSUM = 0.0
337       DO 350 I=1,NUMBG
338 C       ASSUME ISOTROPIC PHOTON EMISSION IN THE LABORATORY SYSTEM
339          CALL GTISO(U1,V1,W1)
340 C       SET THE PHOTON EXIT PARAMETERS
341          UP=U1
342          VP=V1
343          WP=W1
344          AGEP=AGE
345          MTP=MT
346 C re-sample photon energy depending on model used CZ 13.8.92
347          IF(IEP.EQ.2) THEN
348             CALL PARTXS(D(LS2),D(LS2),EOLD,SIGP,EP1)
349             IF(EP1.GT.0.0) EP=EP1
350          ENDIF
351          IF(IEP.EQ.1) THEN
352             CALL SECEGY(EP1,D(LS3),EOLD,D(LS3))
353             IF(EP1.GT.0.0) EP=EP1
354          ENDIF
355          EPSUM = EPSUM+EP
356 C check for energy conservation
357          IF(EPSUM.GT.EPTOT.OR.I.EQ.NUMBG) EP = EPTOT-EPSUM+EP
358 C       STORE THE PHOTON
359          CALL STOPAR(IDGAMA,NGAMA)
360 C end photon production when energy is used up  CZ 13.8.92
361          IF(EPSUM.GT.EPTOT) GOTO 360
362   350 CONTINUE
363   360 RETURN
364   370 WRITE(IOUT,10000)
365 10000 FORMAT(' PHOTON: THE PHOTON PRODUCTION ',
366      +       'CROSS SECTION DATA WAS NOT FOUND (L1=0)')
367       GOTO 390
368   380 WRITE(IOUT,10100)
369 10100 FORMAT(' PHOTON: NO SECONDARY ENERGY ',
370      +   'DISTRIBUTION WAS FOUND FOR THE CONTINUUM REACTION CHOSEN')
371   390 WRITE(6,*) ' CALOR: ERROR in PHOTON ===> STOP '
372       STOP
373       END