]> git.uio.no Git - u/mrichter/AliRoot.git/blame - TFluka/fluscw_1mevn.f
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Correction in ion definition.
[u/mrichter/AliRoot.git] / TFluka / fluscw_1mevn.f
CommitLineData
93d28366 1*$ CREATE FLUSCW.FOR
2*COPY FLUSCW
3* *
4*=== fluscw ===========================================================*
5* *
6 DOUBLE PRECISION FUNCTION FLUSCW
7 #(IJ,PLA,TXX,TYY,TZZ,WEE,XX,YY,ZZ,NREG,IOLREG,LLO,NSURF)
8
9 INCLUDE '(DBLPRC)'
10 INCLUDE '(DIMPAR)'
11 INCLUDE '(IOUNIT)'
12 SAVE
13
14 INCLUDE '(PAPROP)'
15 INCLUDE '(USRBIN)'
16 INCLUDE '(USRBDX)'
17 INCLUDE '(USRTRC)'
18 INCLUDE '(SCOHLP)'
19*
20*----------------------------------------------------------------------*
21* *
22* This functions returns neutron, proton and pion displacement damage *
23* weight factors. *
24* *
25*----------------------------------------------------------------------*
26* *
27* Input variables: *
28* *
29* Ij = (generalized) particle code *
30* Pla = particle momentum (if > 0), or kinetic energy (if <0 )*
31* Txx,yy,zz = particle direction cosines *
32* Wee = particle weight *
33* Xx,Yy,Zz = position *
34* Nreg = (new) region number *
35* Iolreg = (old) region number *
36* Llo = particle generation *
37* Nsurf = transport flag (ignore!) *
38* *
39* Output variables: *
40* *
41* Fluscw = factor the scored amount will be multiplied by *
42* Lsczer = logical flag, if true no amount will be scored *
43* regardless of Fluscw *
44* *
45* Useful variables (common SCOHLP): *
46* *
47* Flux like binnings/estimators (Fluscw): *
48* ISCRNG = 1 --> Boundary crossing estimator *
49* ISCRNG = 2 --> Track length binning *
50* ISCRNG = 3 --> Track length estimator *
51* ISCRNG = 4 --> Collision density estimator *
52* ISCRNG = 5 --> Yield estimator *
53* JSCRNG = # of the binning/estimator *
54* *
55*----------------------------------------------------------------------*
56*
57 LOGICAL LFIRST
58
59
60 DATA LFIRST /.TRUE./
61*
62* the default is unit weight
63 FLUSCW = ONEONE
64 LSCZER = .FALSE.
65*
66* skip all scorings other than tracklength
67 IF (ISCRNG.NE.2) RETURN
68* skip all particles other than protons, pions and neutrons
69 IF ((IJ.NE.1).AND.(IJ.NE.13).AND.(IJ.NE.14).AND.(IJ.NE.8)) RETURN
70*
71* read displacement damage factors
72 IF (LFIRST) THEN
73 WRITE(LUNOUT,1000)
74 1000 FORMAT(1X,'FLUSCW: direct conversion of neutron fluence to',
75 & ' 1 MeV neutron equivalent requested')
76 WRITE(LUNOUT,*) ' neutrons'
77 OPEN(19,FILE='disdam_n.dat',STATUS='UNKNOWN')
78 CALL SKIP(19,7)
79 CALL RDXSC(19,LUNOUT,IDNDAM,-3)
80 CLOSE(19)
81 WRITE(LUNOUT,*) ' protons '
82 OPEN(19,FILE='disdam_p.dat',STATUS='UNKNOWN')
83 CALL SKIP(19,7)
84 CALL RDXSC(19,LUNOUT,IDPDAM,-2)
85 CLOSE(19)
86 WRITE(LUNOUT,*) ' pions '
87 OPEN(19,FILE='disdam_pi.dat',STATUS='UNKNOWN')
88 CALL SKIP(19,7)
89 CALL RDXSC(19,LUNOUT,IDODAM,-2)
90 CLOSE(19)
91 LFIRST = .FALSE.
92 ENDIF
93*
94* should be always called with ekin ( pla < 0 ),
95* but we leave the check for the moment..
96 IF (PLA.LT.0.0D0) THEN
97 EKIN = ABS(PLA)
98 ELSEIF (PLA.GT.0.0D0) THEN
99 EKIN = SQRT(PLA**2+AM(IJ)**2)-AM(IJ)
100 ELSE
101 RETURN
102 ENDIF
103*
104* calculate the weight
105 IF (IJ.EQ.1) THEN
106 FLUSCW = XSECT(IDPDAM,EKIN,0)
107 ELSEIF ((IJ.EQ.13).OR.(IJ.EQ.14)) THEN
108 FLUSCW = XSECT(IDODAM,EKIN,0)
109 ELSEIF (IJ.EQ.8) THEN
110 FLUSCW = XSECT(IDNDAM,EKIN,0)
111 ELSE
112 STOP ' FLUSCW: inconsistent IJ '
113 ENDIF
114
115 RETURN
116 END
117*
118*===rdxsc==============================================================*
119*
120 SUBROUTINE RDXSC(LINP,LCHK,IDXSEC,MODE)
121
122************************************************************************
123* LINP logical input unit *
124* LCHK > 0 cross section data are plotted to unit LCHK *
125* in equidistant log. binning using double-log. *
126* interpolation *
127* (otherwise they are not plotted) *
128* IDXSEC index of stored cross section in common *
129* |MODE| = 1 cross section data given as histogram *
130* i.e. E_lo(i) sigma(i) *
131* E_hi(i) sigma(i) *
132* with increasing energy *
133* = 2 cross section data given for bin averages with *
134* increasing energy *
135* = 3 cross section data given for bin averages with *
136* decreasing energy *
137* if MODE < 0 the energy unit is assumed to be MeV (otherwise GeV) *
138************************************************************************
139
140 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
141 SAVE
142
143 PARAMETER (MAXSDT = 10,
144 & MAXSBI = 1400)
145 COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
146 & NXSBIN(MAXSDT),NXSDT
147 DIMENSION TMPXS(2,MAXSBI)
148
149 LOGICAL LFIRST
150 DATA LFIRST /.TRUE./
151
152 IF (LFIRST) THEN
153 NXSDT = 0
154 DO 1 IDT=1,MAXSDT
155 NXSBIN(IDT) = 0
156 DO 2 IBIN=1,MAXSBI
157 XSEAV(IDT,IBIN) = 0.0D0
158 XS(IDT,IBIN) = 0.0D0
159 2 CONTINUE
160 1 CONTINUE
161 LFIRST = .FALSE.
162 ENDIF
163
164 NXSDT = NXSDT+1
165
166 NEBIN = 0
167 10 CONTINUE
168 NEBIN = NEBIN+1
169 IF (NEBIN.GT.MAXSBI) STOP ' nebin > maxsbi !!!'
170 IF (IABS(MODE).EQ.1) THEN
171 READ(LINP,*,END=9000) ELOW,XS(NXSDT,NEBIN)
172 READ(LINP,*) ELHI,XS(NXSDT,NEBIN)
173 XSEAV(NXSDT,NEBIN) = SQRT(ELOW*ELHI)
174 IF (MODE.LT.0) XSEAV(NXSDT,NEBIN) = 1.D-3*XSEAV(NXSDT,NEBIN)
175 ELSEIF (IABS(MODE).EQ.2) THEN
176 READ(LINP,*,END=9000) XSEAV(NXSDT,NEBIN),XS(NXSDT,NEBIN)
177 IF (MODE.LT.0) XSEAV(NXSDT,NEBIN) = 1.D-3*XSEAV(NXSDT,NEBIN)
178 ELSEIF (IABS(MODE).EQ.3) THEN
179 READ(LINP,*,END=9000) TMPXS(1,NEBIN),TMPXS(2,NEBIN)
180 IF (MODE.LT.0) TMPXS(1,NEBIN) = 1.D-3*TMPXS(1,NEBIN)
181 ELSE
182 STOP ' RDXSC: unsupported mode ! '
183 ENDIF
184 GOTO 10
185 9000 CONTINUE
186
187 IF (IABS(MODE).EQ.3) THEN
188 DO 3 I=1,NEBIN-1
189 IDX = NEBIN-I
190 XSEAV(NXSDT,IDX) = TMPXS(1,I)
191 XS(NXSDT,IDX) = TMPXS(2,I)
192 3 CONTINUE
193 ENDIF
194
195 NXSBIN(NXSDT) = NEBIN-1
196 IDXSEC = NXSDT
197
198 IF (LCHK.GT.0) THEN
199 AELO = LOG10(XSEAV(IDXSEC,1))
200 AEHI = LOG10(XSEAV(IDXSEC,NXSBIN(IDXSEC)))
201 DAE = (AEHI-AELO)/DBLE(NXSBIN(IDXSEC))
202 DO 4 I=1,NXSBIN(IDXSEC)+1
203 AE = AELO+DBLE(I-1)*DAE
204 E = 10.0D0**AE
205 WRITE(LCHK,'(1X,2E15.5)') E,XSECT(IDXSEC,E,0)
206 4 CONTINUE
207 ENDIF
208
209 RETURN
210 END
211*
212*===xsect==============================================================*
213*
214 DOUBLE PRECISION FUNCTION XSECT(IDXSEC,E,MODE)
215
216************************************************************************
217* IDXSEC index of stored cross section in common *
218* E energy *
219************************************************************************
220
221 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
222 SAVE
223
224 PARAMETER (MAXSDT = 10,
225 & MAXSBI = 1400)
226 COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
227 & NXSBIN(MAXSDT),NXSDT
228
229 IF (E.LE.XSEAV(IDXSEC,1)) THEN
230 CS = 0.0D0
231 ELSEIF (E.GT.XSEAV(IDXSEC,NXSBIN(IDXSEC))) THEN
232 N = NXSBIN(IDXSEC)-1
233 CS = XSCINT(IDXSEC,E,N)
234 ELSE
235 DO 1 J=1,NXSBIN(IDXSEC)-1
236 IF ((E.GT.XSEAV(IDXSEC,J)).AND.(E.LE.XSEAV(IDXSEC,J+1)))
237 & THEN
238 CS = XSCINT(IDXSEC,E,J)
239 GOTO 2
240 ENDIF
241 1 CONTINUE
242 STOP ' xsection value not found '
243 2 CONTINUE
244 ENDIF
245 XSECT = CS
246
247 RETURN
248 END
249*
250*===xscint=============================================================*
251*
252 DOUBLE PRECISION FUNCTION XSCINT(IDXSEC,E,J)
253
254 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
255 SAVE
256
257 PARAMETER (MAXSDT = 10,
258 & MAXSBI = 1400)
259 COMMON /CSECT/ XSEAV(MAXSDT,MAXSBI),XS(MAXSDT,MAXSBI),
260 & NXSBIN(MAXSDT),NXSDT
261
262 FAC = (LOG10(E) -LOG10(XSEAV(IDXSEC,J)))/
263 & (LOG10(XSEAV(IDXSEC,J+1))-LOG10(XSEAV(IDXSEC,J)))
264 XSCINT = LOG10(XS(IDXSEC,J))
265 & +(LOG10(XS(IDXSEC,J+1))-LOG10(XS(IDXSEC,J)))*FAC
266 XSCINT = 10**(XSCINT)
267
268 RETURN
269 END
270*
271*===skip===============================================================*
272*
273 SUBROUTINE SKIP(LUNIT,NSKIP)
274
275 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
276 SAVE
277 PARAMETER (LOUT=6,LLOOK=9)
278
279 CHARACTER*132 ACARD
280
281 IF (NSKIP.EQ.0) RETURN
282
283 DO 1 K=1,NSKIP
284 READ(LUNIT,'(A132)') ACARD
285 1 CONTINUE
286
287 RETURN
288 END