1 | // |
---|
2 | // ******************************************************************** |
---|
3 | // * License and Disclaimer * |
---|
4 | // * * |
---|
5 | // * The Geant4 software is copyright of the Copyright Holders of * |
---|
6 | // * the Geant4 Collaboration. It is provided under the terms and * |
---|
7 | // * conditions of the Geant4 Software License, included in the file * |
---|
8 | // * LICENSE and available at http://cern.ch/geant4/license . These * |
---|
9 | // * include a list of copyright holders. * |
---|
10 | // * * |
---|
11 | // * Neither the authors of this software system, nor their employing * |
---|
12 | // * institutes,nor the agencies providing financial support for this * |
---|
13 | // * work make any representation or warranty, express or implied, * |
---|
14 | // * regarding this software system or assume any liability for its * |
---|
15 | // * use. Please see the license in the file LICENSE and URL above * |
---|
16 | // * for the full disclaimer and the limitation of liability. * |
---|
17 | // * * |
---|
18 | // * This code implementation is the result of the scientific and * |
---|
19 | // * technical work of the GEANT4 collaboration. * |
---|
20 | // * By using, copying, modifying or distributing the software (or * |
---|
21 | // * any work based on the software) you agree to acknowledge its * |
---|
22 | // * use in resulting scientific publications, and indicate your * |
---|
23 | // * acceptance of all terms of the Geant4 Software license. * |
---|
24 | // ******************************************************************** |
---|
25 | // |
---|
26 | // |
---|
27 | // $Id: G4QIsotope.cc,v 1.11 2009/02/23 09:49:24 mkossov Exp $ |
---|
28 | // GEANT4 tag $Name: geant4-09-03-beta-cand-01 $ |
---|
29 | // |
---|
30 | // ---------------- G4QIsotope class ---------------- |
---|
31 | // by Mikhail Kossov, December 2003. |
---|
32 | // class G4QIsotope gives Isotopes for Elements (CHIPS branch of G4) |
---|
33 | // ------------------------------------------------------------------ |
---|
34 | // **************************************************************************************** |
---|
35 | // ********** This CLASS is temporary moved from the photolepton_hadron directory ********* |
---|
36 | // ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ****** |
---|
37 | // **************************************************************************************** |
---|
38 | // |
---|
39 | // 1 2 3 4 5 6 7 8 9 |
---|
40 | //34567890123456789012345678901234567890123456789012345678901234567890123456789012345678901 |
---|
41 | // ------------------------------------------------------------------------------ |
---|
42 | // Short description: containes the natural abundance DB for isotops and permits |
---|
43 | // new artificial materials with unnatural abundance (e.g. enreached Deuterium). |
---|
44 | // Uses isotope cross-sections for calculation of the mean cross-sections for the |
---|
45 | // Element (fixed Z). |
---|
46 | // ------------------------------------------------------------------------------ |
---|
47 | |
---|
48 | |
---|
49 | //#define debug |
---|
50 | //#define cdebug |
---|
51 | //#define pdebug |
---|
52 | //#define ppdebug |
---|
53 | //#define sdebug |
---|
54 | |
---|
55 | #include "G4QIsotope.hh" |
---|
56 | #include <cmath> |
---|
57 | using namespace std; |
---|
58 | |
---|
59 | vector<vector<pair<G4int,G4double>*>*>G4QIsotope::natElements;//NaturalElems |
---|
60 | vector<vector<pair<G4int,G4double>*>*>G4QIsotope::natSumAbund;//NaturalSumAb |
---|
61 | vector<vector<pair<G4int,G4double>*>*>G4QIsotope::natIsoCrosS;//CSOfNatElems |
---|
62 | vector<pair<G4int,vector<pair<G4int,G4double>*>*>*>G4QIsotope::newElems; |
---|
63 | vector<pair<G4int,vector<pair<G4int,G4double>*>*>*>G4QIsotope::newSumAb; |
---|
64 | vector<pair<G4int,vector<pair<G4int,G4double>*>*>*>G4QIsotope::newIsoCS; |
---|
65 | |
---|
66 | G4QIsotope::G4QIsotope() |
---|
67 | { |
---|
68 | #ifdef cdebug |
---|
69 | G4cout<<"G4QIsotope::Constructor is called"<<G4endl; |
---|
70 | #endif |
---|
71 | vector<vector<pair<G4int,G4double> >*> natEl; |
---|
72 | #ifdef cdebug |
---|
73 | G4cout<<"G4QIsotope::Constructor natEl is booked"<<G4endl; |
---|
74 | #endif |
---|
75 | vector<pair<G4int,G4double> >*a0=new vector<pair<G4int,G4double> >; |
---|
76 | #ifdef cdebug |
---|
77 | G4cout<<"G4QIsotope::Constructor a0 is booked"<<G4endl; |
---|
78 | #endif |
---|
79 | a0->push_back(make_pair(1,1.)); |
---|
80 | #ifdef cdebug |
---|
81 | G4cout<<"G4QIsotope::Constructor a0 is filled by a pair"<<G4endl; |
---|
82 | #endif |
---|
83 | natEl.push_back(a0); |
---|
84 | #ifdef cdebug |
---|
85 | G4cout<<"G4QIsotope::Constructor a0 is filled in natEl"<<G4endl; |
---|
86 | #endif |
---|
87 | // If an error is found in this initialization, please, correct the simple tree below |
---|
88 | vector<pair<G4int,G4double> >*a1=new vector<pair<G4int,G4double> >; |
---|
89 | a1->push_back(make_pair(0,.99985)); |
---|
90 | a1->push_back(make_pair(1,1.)); |
---|
91 | natEl.push_back(a1); |
---|
92 | vector<pair<G4int,G4double> >*a2=new vector<pair<G4int,G4double> >; |
---|
93 | a2->push_back(make_pair(2,.999999863)); |
---|
94 | a2->push_back(make_pair(1,1.)); |
---|
95 | natEl.push_back(a2); |
---|
96 | vector<pair<G4int,G4double> >*a3=new vector<pair<G4int,G4double> >; |
---|
97 | a3->push_back(make_pair(4,.925)); |
---|
98 | a3->push_back(make_pair(3,1.)); |
---|
99 | natEl.push_back(a3); |
---|
100 | vector<pair<G4int,G4double> >*a4=new vector<pair<G4int,G4double> >; |
---|
101 | a4->push_back(make_pair(5,1.)); |
---|
102 | natEl.push_back(a4); |
---|
103 | vector<pair<G4int,G4double> >*a5=new vector<pair<G4int,G4double> >; |
---|
104 | a5->push_back(make_pair(6,.801)); |
---|
105 | a5->push_back(make_pair(5,1.)); |
---|
106 | natEl.push_back(a5); |
---|
107 | vector<pair<G4int,G4double> >*a6=new vector<pair<G4int,G4double> >; |
---|
108 | a6->push_back(make_pair(6,.989)); |
---|
109 | a6->push_back(make_pair(7,1.)); |
---|
110 | natEl.push_back(a6); |
---|
111 | vector<pair<G4int,G4double> >*a7=new vector<pair<G4int,G4double> >; |
---|
112 | a7->push_back(make_pair(7,.9963)); |
---|
113 | a7->push_back(make_pair(8,1.)); |
---|
114 | natEl.push_back(a7); |
---|
115 | vector<pair<G4int,G4double> >*a8=new vector<pair<G4int,G4double> >; |
---|
116 | a8->push_back(make_pair(8,.9976)); |
---|
117 | a8->push_back(make_pair(10,.9996)); |
---|
118 | a8->push_back(make_pair(9,1.)); |
---|
119 | natEl.push_back(a8); |
---|
120 | vector<pair<G4int,G4double> >*a9=new vector<pair<G4int,G4double> >; |
---|
121 | a9->push_back(make_pair(10,1.)); |
---|
122 | natEl.push_back(a9); |
---|
123 | vector<pair<G4int,G4double> >*b0=new vector<pair<G4int,G4double> >; |
---|
124 | b0->push_back(make_pair(10,.9948)); |
---|
125 | b0->push_back(make_pair(11,.9975)); |
---|
126 | b0->push_back(make_pair(12,1.)); |
---|
127 | natEl.push_back(b0); |
---|
128 | vector<pair<G4int,G4double> >*b1=new vector<pair<G4int,G4double> >; |
---|
129 | b1->push_back(make_pair(12,1.)); |
---|
130 | natEl.push_back(b1); |
---|
131 | vector<pair<G4int,G4double> >*b2=new vector<pair<G4int,G4double> >; |
---|
132 | b2->push_back(make_pair(12,.7899)); |
---|
133 | b2->push_back(make_pair(13,.8899)); |
---|
134 | b2->push_back(make_pair(14,1.)); |
---|
135 | natEl.push_back(b2); |
---|
136 | vector<pair<G4int,G4double> >*b3=new vector<pair<G4int,G4double> >; |
---|
137 | b3->push_back(make_pair(14,1.)); |
---|
138 | natEl.push_back(b3); |
---|
139 | vector<pair<G4int,G4double> >*b4=new vector<pair<G4int,G4double> >; |
---|
140 | b4->push_back(make_pair(14,.9223)); |
---|
141 | b4->push_back(make_pair(15,.969)); |
---|
142 | b4->push_back(make_pair(16,1.)); |
---|
143 | natEl.push_back(b4); |
---|
144 | vector<pair<G4int,G4double> >*b5=new vector<pair<G4int,G4double> >; |
---|
145 | b5->push_back(make_pair(16,1.)); |
---|
146 | natEl.push_back(b5); |
---|
147 | vector<pair<G4int,G4double> >*b6=new vector<pair<G4int,G4double> >; |
---|
148 | b6->push_back(make_pair(16,.9502)); |
---|
149 | b6->push_back(make_pair(18,.9923)); |
---|
150 | b6->push_back(make_pair(17,.9998)); |
---|
151 | b6->push_back(make_pair(20,1.)); |
---|
152 | natEl.push_back(b6); |
---|
153 | vector<pair<G4int,G4double> >*b7=new vector<pair<G4int,G4double> >; |
---|
154 | b7->push_back(make_pair(18,.7577)); |
---|
155 | b7->push_back(make_pair(20,1.)); |
---|
156 | natEl.push_back(b7); |
---|
157 | vector<pair<G4int,G4double> >*b8=new vector<pair<G4int,G4double> >; |
---|
158 | b8->push_back(make_pair(22,.996)); |
---|
159 | b8->push_back(make_pair(18,.99937)); |
---|
160 | b8->push_back(make_pair(20,1.)); |
---|
161 | natEl.push_back(b8); |
---|
162 | vector<pair<G4int,G4double> >*b9=new vector<pair<G4int,G4double> >; |
---|
163 | b9->push_back(make_pair(20,.932581)); |
---|
164 | b9->push_back(make_pair(22,.999883)); |
---|
165 | b9->push_back(make_pair(21,1.)); |
---|
166 | natEl.push_back(b9); |
---|
167 | vector<pair<G4int,G4double> >*c0=new vector<pair<G4int,G4double> >; |
---|
168 | c0->push_back(make_pair(20,.96941)); |
---|
169 | c0->push_back(make_pair(24,.99027)); |
---|
170 | c0->push_back(make_pair(22,.99674)); |
---|
171 | c0->push_back(make_pair(28,.99861)); |
---|
172 | c0->push_back(make_pair(23,.99996)); |
---|
173 | c0->push_back(make_pair(26,1.)); |
---|
174 | natEl.push_back(c0); |
---|
175 | vector<pair<G4int,G4double> >*c1=new vector<pair<G4int,G4double> >; |
---|
176 | c1->push_back(make_pair(24,1.)); |
---|
177 | natEl.push_back(c1); |
---|
178 | vector<pair<G4int,G4double> >*c2=new vector<pair<G4int,G4double> >; |
---|
179 | c2->push_back(make_pair(1,.738)); |
---|
180 | c2->push_back(make_pair(1,.818)); |
---|
181 | c2->push_back(make_pair(1,.891)); |
---|
182 | c2->push_back(make_pair(1,.946)); |
---|
183 | c2->push_back(make_pair(1,1.)); |
---|
184 | natEl.push_back(c2); |
---|
185 | vector<pair<G4int,G4double> >*c3=new vector<pair<G4int,G4double> >; |
---|
186 | c3->push_back(make_pair(28,.9975)); |
---|
187 | c3->push_back(make_pair(27,1.)); |
---|
188 | natEl.push_back(c3); |
---|
189 | vector<pair<G4int,G4double> >*c4=new vector<pair<G4int,G4double> >; |
---|
190 | c4->push_back(make_pair(28,.8379)); |
---|
191 | c4->push_back(make_pair(29,.9329)); |
---|
192 | c4->push_back(make_pair(26,.97635)); |
---|
193 | c4->push_back(make_pair(30,1.)); |
---|
194 | natEl.push_back(c4); |
---|
195 | vector<pair<G4int,G4double> >*c5=new vector<pair<G4int,G4double> >; |
---|
196 | c5->push_back(make_pair(30,1.)); |
---|
197 | natEl.push_back(c5); |
---|
198 | vector<pair<G4int,G4double> >*c6=new vector<pair<G4int,G4double> >; |
---|
199 | c6->push_back(make_pair(30,.9172)); |
---|
200 | c6->push_back(make_pair(28,.9762)); |
---|
201 | c6->push_back(make_pair(31,.9972)); |
---|
202 | c6->push_back(make_pair(32,1.)); |
---|
203 | natEl.push_back(c6); |
---|
204 | vector<pair<G4int,G4double> >*c7=new vector<pair<G4int,G4double> >; |
---|
205 | c7->push_back(make_pair(32,1.)); |
---|
206 | natEl.push_back(c7); |
---|
207 | vector<pair<G4int,G4double> >*c8=new vector<pair<G4int,G4double> >; |
---|
208 | c8->push_back(make_pair(30,.68077)); |
---|
209 | c8->push_back(make_pair(32,.943)); |
---|
210 | c8->push_back(make_pair(34,.97934)); |
---|
211 | c8->push_back(make_pair(33,.99074)); |
---|
212 | c8->push_back(make_pair(36,1.)); |
---|
213 | natEl.push_back(c8); |
---|
214 | vector<pair<G4int,G4double> >*c9=new vector<pair<G4int,G4double> >; |
---|
215 | c9->push_back(make_pair(34,.6917)); |
---|
216 | c9->push_back(make_pair(36,1.)); |
---|
217 | natEl.push_back(c9); |
---|
218 | vector<pair<G4int,G4double> >*d0=new vector<pair<G4int,G4double> >; |
---|
219 | d0->push_back(make_pair(34,.486)); |
---|
220 | d0->push_back(make_pair(36,.765)); |
---|
221 | d0->push_back(make_pair(38,.953)); |
---|
222 | d0->push_back(make_pair(37,.994)); |
---|
223 | d0->push_back(make_pair(40,1.)); |
---|
224 | natEl.push_back(d0); |
---|
225 | vector<pair<G4int,G4double> >*d1=new vector<pair<G4int,G4double> >; |
---|
226 | d1->push_back(make_pair(38,.60108)); |
---|
227 | d1->push_back(make_pair(40,1.)); |
---|
228 | natEl.push_back(d1); |
---|
229 | vector<pair<G4int,G4double> >*d2=new vector<pair<G4int,G4double> >; |
---|
230 | d2->push_back(make_pair(42,.3594)); |
---|
231 | d2->push_back(make_pair(40,.6360)); |
---|
232 | d2->push_back(make_pair(38,.8484)); |
---|
233 | d2->push_back(make_pair(41,.9256)); |
---|
234 | d2->push_back(make_pair(44,1.)); |
---|
235 | natEl.push_back(d2); |
---|
236 | vector<pair<G4int,G4double> >*d3=new vector<pair<G4int,G4double> >; |
---|
237 | d3->push_back(make_pair(42,1.)); |
---|
238 | natEl.push_back(d3); |
---|
239 | vector<pair<G4int,G4double> >*d4=new vector<pair<G4int,G4double> >; |
---|
240 | d4->push_back(make_pair(46,.4961)); |
---|
241 | d4->push_back(make_pair(44,.7378)); |
---|
242 | d4->push_back(make_pair(42,.8274)); |
---|
243 | d4->push_back(make_pair(48,.9148)); |
---|
244 | d4->push_back(make_pair(43,.9911)); |
---|
245 | d4->push_back(make_pair(40,1.)); |
---|
246 | natEl.push_back(d4); |
---|
247 | vector<pair<G4int,G4double> >*d5=new vector<pair<G4int,G4double> >; |
---|
248 | d5->push_back(make_pair(44,.5069)); |
---|
249 | d5->push_back(make_pair(46,1.)); |
---|
250 | natEl.push_back(d5); |
---|
251 | vector<pair<G4int,G4double> >*d6=new vector<pair<G4int,G4double> >; |
---|
252 | d6->push_back(make_pair(48,.57)); |
---|
253 | d6->push_back(make_pair(50,.743)); |
---|
254 | d6->push_back(make_pair(46,.859)); |
---|
255 | d6->push_back(make_pair(47,.974)); |
---|
256 | d6->push_back(make_pair(44,.9965)); |
---|
257 | d6->push_back(make_pair(42,1.)); |
---|
258 | natEl.push_back(d6); |
---|
259 | vector<pair<G4int,G4double> >*d7=new vector<pair<G4int,G4double> >; |
---|
260 | d7->push_back(make_pair(48,.7217)); |
---|
261 | d7->push_back(make_pair(50,1.)); |
---|
262 | natEl.push_back(d7); |
---|
263 | vector<pair<G4int,G4double> >*d8=new vector<pair<G4int,G4double> >; |
---|
264 | d8->push_back(make_pair(50,.8258)); |
---|
265 | d8->push_back(make_pair(48,.9244)); |
---|
266 | d8->push_back(make_pair(49,.9944)); |
---|
267 | d8->push_back(make_pair(46,1.)); |
---|
268 | natEl.push_back(d8); |
---|
269 | vector<pair<G4int,G4double> >*d9=new vector<pair<G4int,G4double> >; |
---|
270 | d9->push_back(make_pair(50,1.)); |
---|
271 | natEl.push_back(d9); |
---|
272 | vector<pair<G4int,G4double> >*e0=new vector<pair<G4int,G4double> >; |
---|
273 | e0->push_back(make_pair(50,.5145)); |
---|
274 | e0->push_back(make_pair(54,.6883)); |
---|
275 | e0->push_back(make_pair(53,.8598)); |
---|
276 | e0->push_back(make_pair(51,.972)); |
---|
277 | e0->push_back(make_pair(56,1.)); |
---|
278 | natEl.push_back(e0); |
---|
279 | vector<pair<G4int,G4double> >*e1=new vector<pair<G4int,G4double> >; |
---|
280 | e1->push_back(make_pair(52,1.)); |
---|
281 | natEl.push_back(e1); |
---|
282 | vector<pair<G4int,G4double> >*e2=new vector<pair<G4int,G4double> >; |
---|
283 | e2->push_back(make_pair(56,.2413)); |
---|
284 | e2->push_back(make_pair(54,.4081)); |
---|
285 | e2->push_back(make_pair(53,.5673)); |
---|
286 | e2->push_back(make_pair(50,.7157)); |
---|
287 | e2->push_back(make_pair(58,.8120)); |
---|
288 | e2->push_back(make_pair(55,.9075)); |
---|
289 | e2->push_back(make_pair(52,1.)); |
---|
290 | natEl.push_back(e2); |
---|
291 | vector<pair<G4int,G4double> >*e3=new vector<pair<G4int,G4double> >; |
---|
292 | e3->push_back(make_pair(55,1.)); |
---|
293 | natEl.push_back(e3); |
---|
294 | vector<pair<G4int,G4double> >*e4=new vector<pair<G4int,G4double> >; |
---|
295 | e4->push_back(make_pair(58,.316)); |
---|
296 | e4->push_back(make_pair(60,.502)); |
---|
297 | e4->push_back(make_pair(57,.673)); |
---|
298 | e4->push_back(make_pair(55,.8)); |
---|
299 | e4->push_back(make_pair(56,.926)); |
---|
300 | e4->push_back(make_pair(52,.9814)); |
---|
301 | e4->push_back(make_pair(54,1.)); |
---|
302 | natEl.push_back(e4); |
---|
303 | vector<pair<G4int,G4double> >*e5=new vector<pair<G4int,G4double> >; |
---|
304 | e5->push_back(make_pair(58,1.)); |
---|
305 | natEl.push_back(e5); |
---|
306 | vector<pair<G4int,G4double> >*e6=new vector<pair<G4int,G4double> >; |
---|
307 | e6->push_back(make_pair(60,.2733)); |
---|
308 | e6->push_back(make_pair(62,.5379)); |
---|
309 | e6->push_back(make_pair(59,.7612)); |
---|
310 | e6->push_back(make_pair(55,.8784)); |
---|
311 | e6->push_back(make_pair(58,.9898)); |
---|
312 | e6->push_back(make_pair(56,1.)); |
---|
313 | natEl.push_back(e6); |
---|
314 | vector<pair<G4int,G4double> >*e7=new vector<pair<G4int,G4double> >; |
---|
315 | e7->push_back(make_pair(60,.51839)); |
---|
316 | e7->push_back(make_pair(62,1.)); |
---|
317 | natEl.push_back(e7); |
---|
318 | vector<pair<G4int,G4double> >*e8=new vector<pair<G4int,G4double> >; |
---|
319 | e8->push_back(make_pair(66,.2873)); |
---|
320 | e8->push_back(make_pair(64,.5286)); |
---|
321 | e8->push_back(make_pair(59,.6566)); |
---|
322 | e8->push_back(make_pair(62,.7815)); |
---|
323 | e8->push_back(make_pair(65,.9037)); |
---|
324 | e8->push_back(make_pair(68,.9786)); |
---|
325 | e8->push_back(make_pair(58,.9911)); |
---|
326 | e8->push_back(make_pair(60,1.)); |
---|
327 | natEl.push_back(e8); |
---|
328 | vector<pair<G4int,G4double> >*e9=new vector<pair<G4int,G4double> >; |
---|
329 | e9->push_back(make_pair(66,.9577)); |
---|
330 | e9->push_back(make_pair(64,1.)); |
---|
331 | natEl.push_back(e9); |
---|
332 | vector<pair<G4int,G4double> >*f0=new vector<pair<G4int,G4double> >; |
---|
333 | f0->push_back(make_pair(70,.3259)); |
---|
334 | f0->push_back(make_pair(68,.5681)); |
---|
335 | f0->push_back(make_pair(66,.7134)); |
---|
336 | f0->push_back(make_pair(69,.7992)); |
---|
337 | f0->push_back(make_pair(67,.8760)); |
---|
338 | f0->push_back(make_pair(74,.9339)); |
---|
339 | f0->push_back(make_pair(72,.9802)); |
---|
340 | f0->push_back(make_pair(62,.9899)); |
---|
341 | f0->push_back(make_pair(64,1.)); |
---|
342 | //f0->push_back(make_pair(64,.9964)); |
---|
343 | //f0->push_back(make_pair(65,1.)); // Nine isotopes is the maximum, so Sn155 is out |
---|
344 | natEl.push_back(f0); |
---|
345 | vector<pair<G4int,G4double> >*f1=new vector<pair<G4int,G4double> >; |
---|
346 | f1->push_back(make_pair(70,.5736)); |
---|
347 | f1->push_back(make_pair(72,1.)); |
---|
348 | natEl.push_back(f1); |
---|
349 | vector<pair<G4int,G4double> >*f2=new vector<pair<G4int,G4double> >; |
---|
350 | f2->push_back(make_pair(78,.3387)); |
---|
351 | f2->push_back(make_pair(76,.6557)); |
---|
352 | f2->push_back(make_pair(74,.8450)); |
---|
353 | f2->push_back(make_pair(73,.9162)); |
---|
354 | f2->push_back(make_pair(72,.9641)); |
---|
355 | f2->push_back(make_pair(70,.9900)); |
---|
356 | f2->push_back(make_pair(71,.99905)); |
---|
357 | f2->push_back(make_pair(68,1.)); |
---|
358 | natEl.push_back(f2); |
---|
359 | vector<pair<G4int,G4double> >*f3=new vector<pair<G4int,G4double> >; |
---|
360 | f3->push_back(make_pair(74,1.)); |
---|
361 | natEl.push_back(f3); |
---|
362 | vector<pair<G4int,G4double> >*f4=new vector<pair<G4int,G4double> >; |
---|
363 | f4->push_back(make_pair(78,.269)); |
---|
364 | f4->push_back(make_pair(75,.533)); |
---|
365 | f4->push_back(make_pair(77,.745)); |
---|
366 | f4->push_back(make_pair(80,.849)); |
---|
367 | f4->push_back(make_pair(82,.938)); |
---|
368 | f4->push_back(make_pair(76,.979)); |
---|
369 | f4->push_back(make_pair(74,.9981)); |
---|
370 | f4->push_back(make_pair(70,.9991)); |
---|
371 | f4->push_back(make_pair(72,1.)); |
---|
372 | natEl.push_back(f4); |
---|
373 | vector<pair<G4int,G4double> >*f5=new vector<pair<G4int,G4double> >; |
---|
374 | f5->push_back(make_pair(78,1.)); |
---|
375 | natEl.push_back(f5); |
---|
376 | vector<pair<G4int,G4double> >*f6=new vector<pair<G4int,G4double> >; |
---|
377 | f6->push_back(make_pair(82,.717)); |
---|
378 | f6->push_back(make_pair(81,.8293)); |
---|
379 | f6->push_back(make_pair(80,.9078)); |
---|
380 | f6->push_back(make_pair(79,.97373)); |
---|
381 | f6->push_back(make_pair(78,.99793)); |
---|
382 | f6->push_back(make_pair(74,.99899)); |
---|
383 | f6->push_back(make_pair(76,1.)); |
---|
384 | natEl.push_back(f6); |
---|
385 | vector<pair<G4int,G4double> >*f7=new vector<pair<G4int,G4double> >; |
---|
386 | f7->push_back(make_pair(82,.999098)); |
---|
387 | f7->push_back(make_pair(81,1.)); |
---|
388 | natEl.push_back(f7); |
---|
389 | vector<pair<G4int,G4double> >*f8=new vector<pair<G4int,G4double> >; |
---|
390 | f8->push_back(make_pair(82,.8843)); |
---|
391 | f8->push_back(make_pair(84,.9956)); |
---|
392 | f8->push_back(make_pair(80,.9981)); |
---|
393 | f8->push_back(make_pair(78,1.)); |
---|
394 | natEl.push_back(f8); |
---|
395 | vector<pair<G4int,G4double> >*f9=new vector<pair<G4int,G4double> >; |
---|
396 | f9->push_back(make_pair(82,1.)); |
---|
397 | natEl.push_back(f9); |
---|
398 | vector<pair<G4int,G4double> >*g0=new vector<pair<G4int,G4double> >; |
---|
399 | g0->push_back(make_pair(82,.2713)); |
---|
400 | g0->push_back(make_pair(84,.5093)); |
---|
401 | g0->push_back(make_pair(86,.6812)); |
---|
402 | g0->push_back(make_pair(83,.8030)); |
---|
403 | g0->push_back(make_pair(85,.8860)); |
---|
404 | g0->push_back(make_pair(88,.9436)); |
---|
405 | g0->push_back(make_pair(90,1.)); |
---|
406 | natEl.push_back(g0); |
---|
407 | vector<pair<G4int,G4double> >*g1=new vector<pair<G4int,G4double> >; |
---|
408 | g1->push_back(make_pair(85,1.)); |
---|
409 | natEl.push_back(g1); |
---|
410 | vector<pair<G4int,G4double> >*g2=new vector<pair<G4int,G4double> >; |
---|
411 | g2->push_back(make_pair(90,.267)); |
---|
412 | g2->push_back(make_pair(92,.494)); |
---|
413 | g2->push_back(make_pair(85,.644)); |
---|
414 | g2->push_back(make_pair(87,.782)); |
---|
415 | g2->push_back(make_pair(86,.895)); |
---|
416 | g2->push_back(make_pair(88,.969)); |
---|
417 | g2->push_back(make_pair(82,1.)); |
---|
418 | natEl.push_back(g2); |
---|
419 | vector<pair<G4int,G4double> >*g3=new vector<pair<G4int,G4double> >; |
---|
420 | g3->push_back(make_pair(90,.522)); |
---|
421 | g3->push_back(make_pair(89,1.)); |
---|
422 | natEl.push_back(g3); |
---|
423 | vector<pair<G4int,G4double> >*g4=new vector<pair<G4int,G4double> >; |
---|
424 | g4->push_back(make_pair(94,.2484)); |
---|
425 | g4->push_back(make_pair(96,.4670)); |
---|
426 | g4->push_back(make_pair(92,.6717)); |
---|
427 | g4->push_back(make_pair(93,.8282)); |
---|
428 | g4->push_back(make_pair(91,.9762)); |
---|
429 | g4->push_back(make_pair(90,.9980)); |
---|
430 | g4->push_back(make_pair(88,1.)); |
---|
431 | natEl.push_back(g4); |
---|
432 | vector<pair<G4int,G4double> >*g5=new vector<pair<G4int,G4double> >; |
---|
433 | g5->push_back(make_pair(94,1.)); |
---|
434 | natEl.push_back(g5); |
---|
435 | vector<pair<G4int,G4double> >*g6=new vector<pair<G4int,G4double> >; |
---|
436 | g6->push_back(make_pair(98,.282)); |
---|
437 | g6->push_back(make_pair(96,.537)); |
---|
438 | g6->push_back(make_pair(97,.786)); |
---|
439 | g6->push_back(make_pair(95,.975)); |
---|
440 | g6->push_back(make_pair(94,.9984)); |
---|
441 | g6->push_back(make_pair(92,.9994)); |
---|
442 | g6->push_back(make_pair(90,1.)); |
---|
443 | natEl.push_back(g6); |
---|
444 | vector<pair<G4int,G4double> >*g7=new vector<pair<G4int,G4double> >; |
---|
445 | g7->push_back(make_pair(98,1.)); |
---|
446 | natEl.push_back(g7); |
---|
447 | vector<pair<G4int,G4double> >*g8=new vector<pair<G4int,G4double> >; |
---|
448 | g8->push_back(make_pair( 98,.3360)); |
---|
449 | g8->push_back(make_pair(100,.6040)); |
---|
450 | g8->push_back(make_pair( 99,.8335)); |
---|
451 | g8->push_back(make_pair(102,.9825)); |
---|
452 | g8->push_back(make_pair( 96,.9986)); |
---|
453 | g8->push_back(make_pair( 94,1.)); |
---|
454 | natEl.push_back(g8); |
---|
455 | vector<pair<G4int,G4double> >*g9=new vector<pair<G4int,G4double> >; |
---|
456 | g9->push_back(make_pair(100,1.)); |
---|
457 | natEl.push_back(g9); |
---|
458 | vector<pair<G4int,G4double> >*h0=new vector<pair<G4int,G4double> >; |
---|
459 | h0->push_back(make_pair(104,.3180)); |
---|
460 | h0->push_back(make_pair(102,.5370)); |
---|
461 | h0->push_back(make_pair(103,.6982)); |
---|
462 | h0->push_back(make_pair(101,.8412)); |
---|
463 | h0->push_back(make_pair(106,.9682)); |
---|
464 | h0->push_back(make_pair(100,.9987)); |
---|
465 | h0->push_back(make_pair( 98,1.)); |
---|
466 | natEl.push_back(h0); |
---|
467 | vector<pair<G4int,G4double> >*h1=new vector<pair<G4int,G4double> >; |
---|
468 | h1->push_back(make_pair(104,.9741)); |
---|
469 | h1->push_back(make_pair(105,1.)); |
---|
470 | natEl.push_back(h1); |
---|
471 | vector<pair<G4int,G4double> >*h2=new vector<pair<G4int,G4double> >; |
---|
472 | h2->push_back(make_pair(108,.35100)); |
---|
473 | h2->push_back(make_pair(106,.62397)); |
---|
474 | h2->push_back(make_pair(105,.81003)); |
---|
475 | h2->push_back(make_pair(107,.94632)); |
---|
476 | h2->push_back(make_pair(104,.99838)); |
---|
477 | h2->push_back(make_pair(102,1.)); |
---|
478 | natEl.push_back(h2); |
---|
479 | vector<pair<G4int,G4double> >*h3=new vector<pair<G4int,G4double> >; |
---|
480 | h3->push_back(make_pair(108,.99988)); |
---|
481 | h3->push_back(make_pair(107,1.)); |
---|
482 | natEl.push_back(h3); |
---|
483 | vector<pair<G4int,G4double> >*h4=new vector<pair<G4int,G4double> >; |
---|
484 | h4->push_back(make_pair(110,.307)); |
---|
485 | h4->push_back(make_pair(112,.593)); |
---|
486 | h4->push_back(make_pair(108,.856)); |
---|
487 | h4->push_back(make_pair(109,.9988)); |
---|
488 | h4->push_back(make_pair(106,1.)); |
---|
489 | natEl.push_back(h4); |
---|
490 | vector<pair<G4int,G4double> >*h5=new vector<pair<G4int,G4double> >; |
---|
491 | h5->push_back(make_pair(112,.626)); |
---|
492 | h5->push_back(make_pair(110,1.)); |
---|
493 | natEl.push_back(h5); |
---|
494 | vector<pair<G4int,G4double> >*h6=new vector<pair<G4int,G4double> >; |
---|
495 | h6->push_back(make_pair(116,.410)); |
---|
496 | h6->push_back(make_pair(114,.674)); |
---|
497 | h6->push_back(make_pair(113,.835)); |
---|
498 | h6->push_back(make_pair(112,.968)); |
---|
499 | h6->push_back(make_pair(111,.984)); |
---|
500 | h6->push_back(make_pair(110,.9998)); |
---|
501 | h6->push_back(make_pair(108,1.)); |
---|
502 | natEl.push_back(h6); |
---|
503 | vector<pair<G4int,G4double> >*h7=new vector<pair<G4int,G4double> >; |
---|
504 | h7->push_back(make_pair(116,.627)); |
---|
505 | h7->push_back(make_pair(114,1.)); |
---|
506 | natEl.push_back(h7); |
---|
507 | vector<pair<G4int,G4double> >*h8=new vector<pair<G4int,G4double> >; |
---|
508 | h8->push_back(make_pair(117,.338)); |
---|
509 | h8->push_back(make_pair(116,.667)); |
---|
510 | h8->push_back(make_pair(118,.920)); |
---|
511 | h8->push_back(make_pair(120,.992)); |
---|
512 | h8->push_back(make_pair(114,.9999)); |
---|
513 | h8->push_back(make_pair(112,1.)); |
---|
514 | natEl.push_back(h8); |
---|
515 | vector<pair<G4int,G4double> >*h9=new vector<pair<G4int,G4double> >; |
---|
516 | h9->push_back(make_pair(118,1.)); |
---|
517 | natEl.push_back(h9); |
---|
518 | vector<pair<G4int,G4double> >*i0=new vector<pair<G4int,G4double> >; |
---|
519 | i0->push_back(make_pair(122,.2986)); |
---|
520 | i0->push_back(make_pair(120,.5296)); |
---|
521 | i0->push_back(make_pair(119,.6983)); |
---|
522 | i0->push_back(make_pair(121,.8301)); |
---|
523 | i0->push_back(make_pair(118,.9298)); |
---|
524 | i0->push_back(make_pair(124,.9985)); |
---|
525 | i0->push_back(make_pair(116,1.)); |
---|
526 | natEl.push_back(i0); |
---|
527 | vector<pair<G4int,G4double> >*i1=new vector<pair<G4int,G4double> >; |
---|
528 | i1->push_back(make_pair(124,.70476)); |
---|
529 | i1->push_back(make_pair(122,1.)); |
---|
530 | natEl.push_back(i1); |
---|
531 | vector<pair<G4int,G4double> >*i2=new vector<pair<G4int,G4double> >; |
---|
532 | i2->push_back(make_pair(126,.524)); |
---|
533 | i2->push_back(make_pair(124,.765)); |
---|
534 | i2->push_back(make_pair(125,.986)); |
---|
535 | i2->push_back(make_pair(122,1.)); |
---|
536 | natEl.push_back(i2); |
---|
537 | vector<pair<G4int,G4double> >*i3=new vector<pair<G4int,G4double> >; |
---|
538 | i3->push_back(make_pair(126,1.)); |
---|
539 | natEl.push_back(i3); |
---|
540 | vector<pair<G4int,G4double> >*i4=new vector<pair<G4int,G4double> >; |
---|
541 | i4->push_back(make_pair(125,1.)); |
---|
542 | natEl.push_back(i4); |
---|
543 | vector<pair<G4int,G4double> >*i5=new vector<pair<G4int,G4double> >; |
---|
544 | i5->push_back(make_pair(136,1.)); |
---|
545 | natEl.push_back(i5); |
---|
546 | vector<pair<G4int,G4double> >*i6=new vector<pair<G4int,G4double> >; |
---|
547 | i6->push_back(make_pair(136,1.)); |
---|
548 | natEl.push_back(i6); |
---|
549 | vector<pair<G4int,G4double> >*i7=new vector<pair<G4int,G4double> >; |
---|
550 | i7->push_back(make_pair(138,1.)); |
---|
551 | natEl.push_back(i7); |
---|
552 | vector<pair<G4int,G4double> >*i8=new vector<pair<G4int,G4double> >; |
---|
553 | i8->push_back(make_pair(138,1.)); |
---|
554 | natEl.push_back(i8); |
---|
555 | vector<pair<G4int,G4double> >*i9=new vector<pair<G4int,G4double> >; |
---|
556 | i9->push_back(make_pair(142,1.)); |
---|
557 | natEl.push_back(i9); |
---|
558 | vector<pair<G4int,G4double> >*j0=new vector<pair<G4int,G4double> >; |
---|
559 | j0->push_back(make_pair(142,1.)); |
---|
560 | natEl.push_back(j0); |
---|
561 | vector<pair<G4int,G4double> >*j1=new vector<pair<G4int,G4double> >; |
---|
562 | j1->push_back(make_pair(140,1.)); |
---|
563 | natEl.push_back(j1); |
---|
564 | vector<pair<G4int,G4double> >*j2=new vector<pair<G4int,G4double> >; |
---|
565 | j2->push_back(make_pair(146,.992745)); |
---|
566 | j2->push_back(make_pair(143,.999945)); |
---|
567 | j2->push_back(make_pair(142,1.)); |
---|
568 | natEl.push_back(j2); |
---|
569 | vector<pair<G4int,G4double> >*j3=new vector<pair<G4int,G4double> >; |
---|
570 | j3->push_back(make_pair(144,1.)); |
---|
571 | natEl.push_back(j3); |
---|
572 | vector<pair<G4int,G4double> >*j4=new vector<pair<G4int,G4double> >; |
---|
573 | j4->push_back(make_pair(150,1.)); |
---|
574 | natEl.push_back(j4); |
---|
575 | vector<pair<G4int,G4double> >*j5=new vector<pair<G4int,G4double> >; |
---|
576 | j5->push_back(make_pair(148,1.)); |
---|
577 | natEl.push_back(j5); |
---|
578 | vector<pair<G4int,G4double> >*j6=new vector<pair<G4int,G4double> >; |
---|
579 | j6->push_back(make_pair(151,1.)); |
---|
580 | natEl.push_back(j6); |
---|
581 | vector<pair<G4int,G4double> >*j7=new vector<pair<G4int,G4double> >; |
---|
582 | j7->push_back(make_pair(150,1.)); |
---|
583 | natEl.push_back(j7); |
---|
584 | vector<pair<G4int,G4double> >*j8=new vector<pair<G4int,G4double> >; |
---|
585 | j8->push_back(make_pair(153,1.)); |
---|
586 | natEl.push_back(j8); |
---|
587 | vector<pair<G4int,G4double> >*j9=new vector<pair<G4int,G4double> >; |
---|
588 | j9->push_back(make_pair(157,1.)); |
---|
589 | natEl.push_back(j9); |
---|
590 | vector<pair<G4int,G4double> >*k0=new vector<pair<G4int,G4double> >; |
---|
591 | k0->push_back(make_pair(157,1.)); |
---|
592 | natEl.push_back(k0); |
---|
593 | vector<pair<G4int,G4double> >*k1=new vector<pair<G4int,G4double> >; |
---|
594 | k1->push_back(make_pair(157,1.)); |
---|
595 | natEl.push_back(k1); |
---|
596 | vector<pair<G4int,G4double> >*k2=new vector<pair<G4int,G4double> >; |
---|
597 | k2->push_back(make_pair(157,1.)); |
---|
598 | natEl.push_back(k2); |
---|
599 | vector<pair<G4int,G4double> >*k3=new vector<pair<G4int,G4double> >; |
---|
600 | k3->push_back(make_pair(157,1.)); |
---|
601 | natEl.push_back(k3); |
---|
602 | vector<pair<G4int,G4double> >*k4=new vector<pair<G4int,G4double> >; |
---|
603 | k4->push_back(make_pair(157,1.)); |
---|
604 | natEl.push_back(k4); |
---|
605 | vector<pair<G4int,G4double> >*k5=new vector<pair<G4int,G4double> >; |
---|
606 | k5->push_back(make_pair(157,1.)); |
---|
607 | natEl.push_back(k5); |
---|
608 | vector<pair<G4int,G4double> >*k6=new vector<pair<G4int,G4double> >; |
---|
609 | k6->push_back(make_pair(157,1.)); |
---|
610 | natEl.push_back(k6); |
---|
611 | vector<pair<G4int,G4double> >*k7=new vector<pair<G4int,G4double> >; |
---|
612 | k7->push_back(make_pair(155,1.)); |
---|
613 | natEl.push_back(k7); |
---|
614 | vector<pair<G4int,G4double> >*k8=new vector<pair<G4int,G4double> >; |
---|
615 | k8->push_back(make_pair(157,1.)); |
---|
616 | natEl.push_back(k8); |
---|
617 | vector<pair<G4int,G4double> >*k9=new vector<pair<G4int,G4double> >; |
---|
618 | k9->push_back(make_pair(157,1.)); |
---|
619 | natEl.push_back(k9); |
---|
620 | // Now fill natElements and natIsoCrossS |
---|
621 | G4int nona=natEl.size(); |
---|
622 | #ifdef cdebug |
---|
623 | G4cout<<"G4QIsotope::Constructor natEl filling is finished nE="<<nona<<G4endl; |
---|
624 | #endif |
---|
625 | for(G4int i=0; i<nona; i++) |
---|
626 | { |
---|
627 | vector<pair<G4int,G4double> >* is=natEl[i]; // Pointer to theElement |
---|
628 | G4int n=is->size(); |
---|
629 | #ifdef cdebug |
---|
630 | G4cout<<"G4QIsotope::Constructor: Element # "<<i<<", nOfIsotopes="<<n<<G4endl; |
---|
631 | #endif |
---|
632 | vector<pair<G4int,G4double>*>*a=new vector<pair<G4int,G4double>*>; |
---|
633 | vector<pair<G4int,G4double>*>*s=new vector<pair<G4int,G4double>*>; |
---|
634 | G4double last=0.; |
---|
635 | if(n) for(G4int j=0; j<n; j++) |
---|
636 | { |
---|
637 | G4int nn =(*is)[j].first; // #ofNeutrons in the isotope |
---|
638 | G4double cur=(*is)[j].second;// value of the summed abundancy |
---|
639 | pair<G4int,G4double>* aP = new pair<G4int,G4double>(nn,cur-last); |
---|
640 | last=cur; // Update the summed value |
---|
641 | pair<G4int,G4double>* sP = new pair<G4int,G4double>((*is)[j]); |
---|
642 | a->push_back(aP); |
---|
643 | s->push_back(sP); |
---|
644 | #ifdef cdebug |
---|
645 | G4cout<<"G4QIsotope::Constructor:Element# "<<i<<", Pair # "<<j<<" is filled"<<G4endl; |
---|
646 | #endif |
---|
647 | } |
---|
648 | natElements.push_back(a); // Fill abundancies for the particular isotope |
---|
649 | natSumAbund.push_back(s); // Fill summes abundancies up to this isotope |
---|
650 | #ifdef cdebug |
---|
651 | G4cout<<"G4QIsotope::Constructor: natElements is filled"<<G4endl; |
---|
652 | #endif |
---|
653 | vector<pair<G4int,G4double>*>*c=new vector<pair<G4int,G4double>*>; |
---|
654 | if(n) for(G4int j=0; j<n; j++) // Cross sections are 0. by default |
---|
655 | { |
---|
656 | pair<G4int,G4double>* cP = new pair<G4int,G4double>((*is)[j].first,0.); |
---|
657 | c->push_back(cP); |
---|
658 | #ifdef cdebug |
---|
659 | G4cout<<"G4QIsotope::Constructor:CrosSecPair i="<<i<<", j="<<j<<" is filled"<<G4endl; |
---|
660 | #endif |
---|
661 | } |
---|
662 | natIsoCrosS.push_back(c); // FillPrototypeCrossSec's (0) for theParticularIsotope |
---|
663 | #ifdef cdebug |
---|
664 | G4cout<<"G4QIsotope::Constructor: natIsoCrosS is filled"<<G4endl; |
---|
665 | #endif |
---|
666 | delete is; |
---|
667 | } |
---|
668 | #ifdef cdebug |
---|
669 | G4cout<<"G4QIsotope::Constructor: is finished"<<G4endl; |
---|
670 | #endif |
---|
671 | } |
---|
672 | |
---|
673 | G4QIsotope::~G4QIsotope() // The QIsotopes are destructed only in theEnd of theJob |
---|
674 | { |
---|
675 | #ifdef debug |
---|
676 | G4cout<<"G4QIsotope::Destructor is called"<<G4endl; |
---|
677 | #endif |
---|
678 | G4int uP=natElements.size(); |
---|
679 | if(uP) for(G4int i=0; i<uP; i++) |
---|
680 | { |
---|
681 | vector<pair<G4int,G4double>*>* curA=natElements[i]; |
---|
682 | G4int nn=curA->size(); // Can not be 0 by definition |
---|
683 | if(nn) for(G4int n=0; n<nn; n++) delete (*curA)[n]; // Delete pair(N,Ab) |
---|
684 | delete curA; // Delet abundancy vector |
---|
685 | vector<pair<G4int,G4double>*>* curS=natSumAbund[i]; |
---|
686 | G4int ns=curS->size(); // Can not be 0 by definition |
---|
687 | if(ns) for(G4int m=0; m<ns; m++) delete (*curS)[m]; // Delete pair(N,Ab) |
---|
688 | delete curS; // Delet abundancy vector |
---|
689 | vector<pair<G4int,G4double>*>* curC=natIsoCrosS[i]; |
---|
690 | G4int nc=curC->size(); // Can not be 0 by definition |
---|
691 | if(nc) for(G4int k=0; k<nc; k++) delete (*curC)[k]; // Delete pair(N,CS) |
---|
692 | delete curC; // Delete cross section vector |
---|
693 | } |
---|
694 | G4int nP=newElems.size(); |
---|
695 | if(nP) for(G4int j=0; j<nP; j++) // LOOP over new UserDefinedElements |
---|
696 | { |
---|
697 | pair<G4int, vector<pair<G4int,G4double>*>* >* nEl= newElems[j]; |
---|
698 | G4int nEn=nEl->second->size(); |
---|
699 | if(nEn) for(G4int k=0; k<nEn; k++) |
---|
700 | { |
---|
701 | pair<G4int,G4double>* curA=(*(nEl->second))[k]; |
---|
702 | delete curA; // Delete vect<pair(N,Abundancy)*> |
---|
703 | } |
---|
704 | delete nEl; // Delete vect<IndZ,vect<pair(N,Ab)*>*> newElementVector |
---|
705 | // |
---|
706 | pair<G4int, vector<pair<G4int,G4double>*>* >* nSA= newSumAb[j]; |
---|
707 | G4int nSn=nSA->second->size(); |
---|
708 | if(nSn) for(G4int n=0; n<nSn; n++) |
---|
709 | { |
---|
710 | pair<G4int,G4double>* curS=(*(nSA->second))[n]; |
---|
711 | delete curS; // Delete vect<pair(N,SumAbund)*> |
---|
712 | } |
---|
713 | delete nSA; // Delete vect<IndZ,vect<pair(N,SA)*>*> newSumAbunVector |
---|
714 | // |
---|
715 | pair<G4int, vector<pair<G4int,G4double>*>* >* nCS= newIsoCS[j]; |
---|
716 | G4int nCn=nCS->second->size(); |
---|
717 | if(nCn) for(G4int m=0; m<nCn; m++) |
---|
718 | { |
---|
719 | pair<G4int,G4double>* curC = (*(nCS->second))[m]; |
---|
720 | delete curC; // Delete vect<pair(N,CrossSect)*> |
---|
721 | } |
---|
722 | delete nCS; // Delete vect<IndZ,vect<pair(N,CS)*>*> newIsoCroSVector |
---|
723 | // |
---|
724 | if(nEn!=nCn) G4cerr<<"*G4QIsotope-WORNING-:#El="<<j<<":nE="<<nEn<<"!=nC="<<nCn<<G4endl; |
---|
725 | if(nEn!=nSn) G4cerr<<"*G4QIsotope-WORNING-:#El="<<j<<":nE="<<nEn<<"!=nS="<<nSn<<G4endl; |
---|
726 | } |
---|
727 | } |
---|
728 | |
---|
729 | // Returns Pointer to the G4QIsotope singletone |
---|
730 | G4QIsotope* G4QIsotope::Get() |
---|
731 | // ================= |
---|
732 | { |
---|
733 | #ifdef pdebug |
---|
734 | G4cout<<"G4QIsotope::Get is called"<<G4endl; |
---|
735 | #endif |
---|
736 | static G4QIsotope theIsotopes; // *** Static body of the G4QIsotope class *** |
---|
737 | return &theIsotopes; |
---|
738 | } |
---|
739 | |
---|
740 | // #ofProtons in stable isotopes with fixed A=Z+N. Returns length and fils VectOfIsotopes |
---|
741 | G4int G4QIsotope::GetProtons(G4int A, vector<G4int>& isoV) |
---|
742 | // ========================================================= |
---|
743 | { |
---|
744 | const G4int nAZ=270; // Dimension of the table |
---|
745 | // Best Z for the given A |
---|
746 | const G4int bestZ[nAZ] = { |
---|
747 | 0, 1, 1, 2, 2, 0, 3, 3, 4, 4, //0 |
---|
748 | 5, 5, 6, 6, 7, 7, 8, 8, 8, 9, //10 |
---|
749 | 10, 10, 10, 11, 12, 12, 12, 13, 14, 14, //20 |
---|
750 | 14, 15, 16, 16, 16, 17, 18, 17, 18, 19, //30 |
---|
751 | 18, 19, 20, 20, 20, 21, 22, 22, 23, 23, //40 |
---|
752 | 22, 23, 24, 24, 26, 25, 26, 26, 28, 27, //50 |
---|
753 | 28, 28, 28, 29, 30, 29, 30, 30, 30, 31, //60 |
---|
754 | 32, 31, 32, 32, 32, 33, 34, 34, 34, 35, //70 |
---|
755 | 34, 35, 36, 36, 36, 37, 39, 36, 38, 39, //80 |
---|
756 | 40, 40, 41, 40, 40, 42, 42, 42, 42, 44, //90 |
---|
757 | 44, 44, 44, 45, 44, 46, 46, 47, 46, 47, //100 |
---|
758 | 48, 48, 48, 48, 48, 49, 50, 50, 50, 50, //110 |
---|
759 | 50, 51, 50, 51, 50, 52, 52, 53, 52, 54, //120 |
---|
760 | 52, 54, 54, 55, 54, 56, 54, 56, 56, 57, //130 |
---|
761 | 58, 59, 60, 60, 60, 60, 60, 62, 62, 62, //140 |
---|
762 | 62, 63, 62, 63, 62, 64, 64, 64, 64, 65, //150 |
---|
763 | 64, 66, 66, 66, 66, 67, 68, 68, 68, 69, //160 |
---|
764 | 68, 70, 70, 70, 70, 71, 70, 72, 72, 72, //170 |
---|
765 | 72, 73, 74, 74, 74, 75, 74, 75, 76, 76, //180 |
---|
766 | 76, 77, 76, 77, 78, 78, 78, 79, 80, 80, //190 |
---|
767 | 80, 80, 80, 81, 80, 81, 82, 82, 82, 83, //200 |
---|
768 | 82, 0, 82, 0, 82, 0, 84, 0, 0, 0, //210 |
---|
769 | 86, 0, 86, 87, 88, 0, 88, 89, 88, 89, //220 |
---|
770 | 89, 91, 90, 0, 92, 92, 0, 93, 92, 94, //230 |
---|
771 | 0, 0, 0, 95, 94, 0, 0, 96, 0, 0, //240 |
---|
772 | 0, 98, 99, 0, 0, 0, 0,100,101,102, //250 |
---|
773 | 103,104,105,106, 0,108,109, 0, 0, 0}; //260 |
---|
774 | // 0 1 2 3 4 5 6 7 8 9 |
---|
775 | // Second candidate |
---|
776 | const G4int secoZ[nAZ] = { |
---|
777 | 0, 0, 0, 1, 0, 0, 0, 4, 0, 0, //0 |
---|
778 | 4, 6, 5, 7, 0, 8, 0, 0, 0, 8, //10 |
---|
779 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //20 |
---|
780 | 0, 0, 15, 0, 0, 16, 0, 0, 0, 0, //30 |
---|
781 | 20, 20, 0, 0, 0, 0, 20, 0, 0, 0, //40 |
---|
782 | 24, 0, 0, 0, 24, 0, 0, 0, 26, 0, //50 |
---|
783 | 27, 0, 0, 0, 28, 0, 0, 0, 0, 0, //60 |
---|
784 | 30, 0, 0, 0, 34, 0, 32, 0, 0, 0, //70 |
---|
785 | 36, 0, 34, 0, 38, 0, 38, 38, 0, 0, //80 |
---|
786 | 0, 0, 42, 0, 42, 0, 44, 0, 44, 0, //90 |
---|
787 | 42, 0, 46, 0, 46, 0, 48, 0, 48, 0, //100 |
---|
788 | 46, 0, 50, 49, 50, 50, 48, 0, 0, 0, //110 |
---|
789 | 52, 0, 52, 0, 52, 0, 54, 0, 54, 0, //120 |
---|
790 | 54, 0, 56, 0, 56, 0, 56, 0, 58, 0, //130 |
---|
791 | 54, 0, 58, 0, 62, 61, 0, 0, 60, 0, //140 |
---|
792 | 60, 0, 64, 0, 64, 0, 66, 0, 66, 0, //150 |
---|
793 | 66, 0, 68, 0, 68, 0, 0, 0, 70, 0, //160 |
---|
794 | 70, 0, 0, 0, 72, 0, 72, 0, 0, 0, //170 |
---|
795 | 74, 0, 0, 0, 76, 0, 76, 76, 0, 0, //180 |
---|
796 | 78, 0, 78, 0, 0, 0, 80, 0, 78, 0, //190 |
---|
797 | 0, 0, 0, 0, 82, 0, 0, 0, 0, 84, //200 |
---|
798 | 84, 0, 83, 0, 83, 0, 0, 0, 0, 0, //210 |
---|
799 | 0, 0, 0, 0, 0, 0, 0, 0, 89, 0, //220 |
---|
800 | 0, 0, 0, 0, 93, 0, 0, 0, 93, 0, //230 |
---|
801 | 0, 0, 0, 0, 0, 0, 0, 97, 0, 0, //240 |
---|
802 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //250 |
---|
803 | 0, 0,107, 0, 0, 0, 0, 0, 0, 0}; //260 |
---|
804 | // 0 1 2 3 4 5 6 7 8 9 |
---|
805 | // Third candidate |
---|
806 | const G4int thrdZ[nAZ] = { |
---|
807 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //0 |
---|
808 | 0, 0, 7, 0, 0, 0, 0, 0, 0, 0, //10 |
---|
809 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //20 |
---|
810 | 0, 0, 0, 0, 0,20, 0, 0, 0, 0, //30 |
---|
811 | 19, 0, 0, 0, 0, 0, 0, 0, 0, 0, //40 |
---|
812 | 23, 0, 0, 0, 0, 0, 0, 0, 0, 0, //50 |
---|
813 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //60 |
---|
814 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //70 |
---|
815 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //80 |
---|
816 | 0, 0,36, 0,38, 0,40, 0, 0, 0, //90 |
---|
817 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //100 |
---|
818 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //110 |
---|
819 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //120 |
---|
820 | 56, 0,50, 0, 0, 0,58, 0,57, 0, //130 |
---|
821 | 0, 0, 0, 0, 0, 0, 0, 0,65, 0, //140 |
---|
822 | 0, 0,66, 0, 0, 0, 0, 0, 0, 0, //150 |
---|
823 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //160 |
---|
824 | 0, 0, 0, 0, 0, 0,71, 0, 0, 0, //170 |
---|
825 | 73, 0, 0, 0, 0, 0, 0, 0, 0, 0, //180 |
---|
826 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //190 |
---|
827 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //200 |
---|
828 | 83, 0,84, 0,84, 0, 0, 0, 0, 0, //210 |
---|
829 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //220 |
---|
830 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //230 |
---|
831 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //240 |
---|
832 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //250 |
---|
833 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; //260 |
---|
834 | //0 1 2 3 4 5 6 7 8 9 |
---|
835 | // Fourth candidate (only two isotopes) |
---|
836 | const G4int quadZ[nAZ] = { |
---|
837 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //0 |
---|
838 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //10 |
---|
839 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //20 |
---|
840 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //30 |
---|
841 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //40 |
---|
842 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //50 |
---|
843 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //60 |
---|
844 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //70 |
---|
845 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //80 |
---|
846 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //90 |
---|
847 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //100 |
---|
848 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //110 |
---|
849 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //120 |
---|
850 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //130 |
---|
851 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //140 |
---|
852 | 0, 0,67, 0, 0, 0, 0, 0, 0, 0, //150 |
---|
853 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //160 |
---|
854 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //170 |
---|
855 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //180 |
---|
856 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //190 |
---|
857 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //200 |
---|
858 | 85, 0, 0, 0, 0, 0, 0, 0, 0, 0, //210 |
---|
859 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //220 |
---|
860 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //230 |
---|
861 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //240 |
---|
862 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, //250 |
---|
863 | 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; //260 |
---|
864 | //0 1 2 3 4 5 6 7 8 9 |
---|
865 | isoV.clear(); // Always cleans up before filling |
---|
866 | if(A>=nAZ) return 0; |
---|
867 | G4int fZ=bestZ[A]; |
---|
868 | if(fZ) |
---|
869 | { |
---|
870 | isoV.push_back(fZ); |
---|
871 | G4int sZ=secoZ[A]; |
---|
872 | if(sZ) |
---|
873 | { |
---|
874 | isoV.push_back(sZ); |
---|
875 | G4int tZ=thrdZ[A]; |
---|
876 | if(tZ) |
---|
877 | { |
---|
878 | isoV.push_back(tZ); |
---|
879 | G4int qZ=quadZ[A]; |
---|
880 | if(qZ) |
---|
881 | { |
---|
882 | isoV.push_back(qZ); |
---|
883 | return 4; |
---|
884 | } |
---|
885 | else return 3; |
---|
886 | } |
---|
887 | else return 2; |
---|
888 | } |
---|
889 | else return 1; |
---|
890 | } |
---|
891 | else return 0; |
---|
892 | } |
---|
893 | |
---|
894 | // Randomize a#ofNeutrons in the Element (independent function @@ can be used for initial) |
---|
895 | G4int G4QIsotope::RandomizeNeutrons(G4int i) |
---|
896 | // ====================================== |
---|
897 | { |
---|
898 | static const G4int nElements = 110; // Max=Meitnerium(Mt)Z=99(starts with Z=0 - neuteron) |
---|
899 | // If an error is found in this simple tree, please, correct the initialization above |
---|
900 | static G4int dN[nElements]= |
---|
901 | { |
---|
902 | // n Be F Al P |
---|
903 | 1, -1, -1, -1, 5, -1, -1, -1, -1, 10, -1, 12, -1, 14, -1, 16, -1, -1, -1, -1, |
---|
904 | // Sc Mn Co As Y |
---|
905 | -1, 24, -1, -1, -1, 30, -1, 32, -1, -1, -1, -1, -1, 42, -1, -1, -1, -1, -1, 50, |
---|
906 | // Nb Tc Rh I Cs Pr |
---|
907 | -1, 52, -1, 55, -1, 58, -1, -1, -1, -1, -1, -1, -1, 74, -1, 78, -1, -1, -1, 82, |
---|
908 | // Pm Tb Ho Tm Au |
---|
909 | -1,-85, -1, -1, -1, 94, -1, 98, -1,100, -1, -1, -1, -1, -1, -1, -1, -1, -1,118, |
---|
910 | // Po At Rn Fr Ra Ac |
---|
911 | -1, -1, -1, 126,-125,-136,-136,-138,-138,-142, |
---|
912 | // Th Pa U Np Pu Am Cm Bk Cf Es |
---|
913 | 142,-140, -1,-144,-150,-148,-151,-150,-153,-153 |
---|
914 | // Fm Md No Lr Rf Db Sg Bh Hs Mt |
---|
915 | -157,-157,-157,-157,-157,-157,-157,-155,-157,-157 |
---|
916 | }; |
---|
917 | #ifdef debug |
---|
918 | G4cout<<"G4QIsotope::RandomizeNeutrons is called Z="<<i<<G4endl; |
---|
919 | #endif |
---|
920 | G4int N=dN[i]; |
---|
921 | if(N==-1) |
---|
922 | { |
---|
923 | G4double rnd=G4UniformRand(); |
---|
924 | if (i<44) // ====== H - Mo |
---|
925 | { |
---|
926 | if (i<23) // ------ H - Ti |
---|
927 | { |
---|
928 | if (i<12) // ______ H - Ne |
---|
929 | { |
---|
930 | if (i<6) // ...... H - B |
---|
931 | { |
---|
932 | if (i<3) |
---|
933 | { |
---|
934 | if(i==1) // H |
---|
935 | { |
---|
936 | if(rnd>.00015) N=0; |
---|
937 | else N=1; |
---|
938 | } |
---|
939 | else // He (2) |
---|
940 | { |
---|
941 | if(rnd>1.37e-6) N=2; |
---|
942 | else N=1; |
---|
943 | } |
---|
944 | } |
---|
945 | else |
---|
946 | { |
---|
947 | if(i==3) // Li |
---|
948 | { |
---|
949 | if(rnd>.075) N=4; |
---|
950 | else N=3; |
---|
951 | } |
---|
952 | else // B (5) |
---|
953 | { |
---|
954 | if(rnd>.199) N=6; |
---|
955 | else N=5; |
---|
956 | } |
---|
957 | } |
---|
958 | } |
---|
959 | else // ...... C - Ne |
---|
960 | { |
---|
961 | if (i<8) |
---|
962 | { |
---|
963 | if(i==6) // C |
---|
964 | { |
---|
965 | if(rnd>.011) N=6; |
---|
966 | else N=7; |
---|
967 | } |
---|
968 | else // N (7) |
---|
969 | { |
---|
970 | if(rnd>.0037) N=7; |
---|
971 | else N=8; |
---|
972 | } |
---|
973 | } |
---|
974 | else |
---|
975 | { |
---|
976 | if(i==8) // O |
---|
977 | { |
---|
978 | if (rnd<.9976) N=8; |
---|
979 | else if(rnd<.9996) N=10; |
---|
980 | else N=9; |
---|
981 | } |
---|
982 | else // Ne (10) |
---|
983 | { |
---|
984 | if (rnd<.9948) N=10; |
---|
985 | else if(rnd<.9975) N=11; |
---|
986 | else N=12; |
---|
987 | } |
---|
988 | } |
---|
989 | } |
---|
990 | } |
---|
991 | else // ______ Mg - Ti |
---|
992 | { |
---|
993 | if (i<18) // ...... Mg - Cl |
---|
994 | { |
---|
995 | if (i<16) |
---|
996 | { |
---|
997 | if(i==12) // Mg |
---|
998 | { |
---|
999 | if (rnd<.7899) N=12; |
---|
1000 | else if(rnd<.8899) N=13; |
---|
1001 | else N=14; |
---|
1002 | } |
---|
1003 | else // Si (14) |
---|
1004 | { |
---|
1005 | if (rnd<.9223) N=14; |
---|
1006 | else if(rnd<.969) N=15; |
---|
1007 | else N=16; |
---|
1008 | } |
---|
1009 | } |
---|
1010 | else |
---|
1011 | { |
---|
1012 | if(i==16) // S |
---|
1013 | { |
---|
1014 | if (rnd<.9502) N=16; |
---|
1015 | else if(rnd<.9923) N=18; |
---|
1016 | else if(rnd<.9998) N=17; |
---|
1017 | else N=20; |
---|
1018 | } |
---|
1019 | else // Cl (17) |
---|
1020 | { |
---|
1021 | if (rnd>.7577) N=18; |
---|
1022 | else N=20; |
---|
1023 | } |
---|
1024 | } |
---|
1025 | } |
---|
1026 | else // ...... Ar - Ti |
---|
1027 | { |
---|
1028 | if (i<20) |
---|
1029 | { |
---|
1030 | if(i==18) // Ar |
---|
1031 | { |
---|
1032 | if (rnd<.996) N=22; |
---|
1033 | else if(rnd<.99937) N=18; |
---|
1034 | else N=20; |
---|
1035 | } |
---|
1036 | else // K (19) |
---|
1037 | { |
---|
1038 | if (rnd<.932581) N=20; |
---|
1039 | else if(rnd<.999883) N=22; |
---|
1040 | else N=21; |
---|
1041 | } |
---|
1042 | } |
---|
1043 | else |
---|
1044 | { |
---|
1045 | if(i==20) // Ca |
---|
1046 | { |
---|
1047 | if (rnd<.96941) N=20; |
---|
1048 | else if(rnd<.99027) N=24; |
---|
1049 | else if(rnd<.99674) N=22; |
---|
1050 | else if(rnd<.99861) N=28; |
---|
1051 | else if(rnd<.99996) N=23; |
---|
1052 | else N=26; |
---|
1053 | } |
---|
1054 | else // Ti (22) |
---|
1055 | { |
---|
1056 | if (rnd<.738) N=26; |
---|
1057 | else if(rnd<.818) N=24; |
---|
1058 | else if(rnd<.891) N=25; |
---|
1059 | else if(rnd<.946) N=27; |
---|
1060 | else N=28; |
---|
1061 | } |
---|
1062 | } |
---|
1063 | } |
---|
1064 | } |
---|
1065 | } |
---|
1066 | else // ------ V - Mo |
---|
1067 | { |
---|
1068 | if (i<32) // ______ V - Ga |
---|
1069 | { |
---|
1070 | if (i<28) // ...... H - Fe |
---|
1071 | { |
---|
1072 | if (i<26) |
---|
1073 | { |
---|
1074 | if(i==23) // V |
---|
1075 | { |
---|
1076 | if (rnd<.9975) N=28; |
---|
1077 | else N=27; |
---|
1078 | } |
---|
1079 | else // Cr (24) |
---|
1080 | { |
---|
1081 | if (rnd<.8379) N=28; |
---|
1082 | else if(rnd<.9329) N=29; |
---|
1083 | else if(rnd<.97635) N=26; |
---|
1084 | else N=30; |
---|
1085 | } |
---|
1086 | } |
---|
1087 | else // Fe (26) |
---|
1088 | { |
---|
1089 | if (rnd<.9172) N=30; |
---|
1090 | else if(rnd<.9762) N=28; |
---|
1091 | else if(rnd<.9972) N=31; |
---|
1092 | else N=32; |
---|
1093 | } |
---|
1094 | } |
---|
1095 | else // ...... Ni - Ga |
---|
1096 | { |
---|
1097 | if (i<30) |
---|
1098 | { |
---|
1099 | if(i==28) // Ni |
---|
1100 | { |
---|
1101 | if (rnd<.68077) N=30; |
---|
1102 | else if(rnd<.943) N=32; |
---|
1103 | else if(rnd<.97934) N=34; |
---|
1104 | else if(rnd<.99074) N=33; |
---|
1105 | else N=36; |
---|
1106 | } |
---|
1107 | else // Cu (29) |
---|
1108 | { |
---|
1109 | if (rnd<.6917) N=34; |
---|
1110 | else N=36; |
---|
1111 | } |
---|
1112 | } |
---|
1113 | else |
---|
1114 | { |
---|
1115 | if(i==30) // Zn |
---|
1116 | { |
---|
1117 | if (rnd<.486) N=34; |
---|
1118 | else if(rnd<.765) N=36; |
---|
1119 | else if(rnd<.953) N=38; |
---|
1120 | else if(rnd<.994) N=37; |
---|
1121 | else N=40; |
---|
1122 | } |
---|
1123 | else // Ga (31) |
---|
1124 | { |
---|
1125 | if (rnd<.60108) N=38; |
---|
1126 | else N=40; |
---|
1127 | } |
---|
1128 | } |
---|
1129 | } |
---|
1130 | } |
---|
1131 | else // ______ Ge - Mo |
---|
1132 | { |
---|
1133 | if (i<37) // ...... H - B |
---|
1134 | { |
---|
1135 | if (i<35) |
---|
1136 | { |
---|
1137 | if(i==32) // Ge |
---|
1138 | { |
---|
1139 | if (rnd<.3594) N=42; |
---|
1140 | else if(rnd<.6360) N=40; |
---|
1141 | else if(rnd<.8484) N=38; |
---|
1142 | else if(rnd<.9256) N=41; |
---|
1143 | else N=44; |
---|
1144 | } |
---|
1145 | else // Se (34) |
---|
1146 | { |
---|
1147 | if (rnd>.4961) N=46; |
---|
1148 | else if(rnd<.7378) N=44; |
---|
1149 | else if(rnd<.8274) N=42; |
---|
1150 | else if(rnd<.9148) N=48; |
---|
1151 | else if(rnd<.9911) N=43; |
---|
1152 | else N=40; |
---|
1153 | } |
---|
1154 | } |
---|
1155 | else |
---|
1156 | { |
---|
1157 | if(i==35) // Br |
---|
1158 | { |
---|
1159 | if (rnd<.5069) N=44; |
---|
1160 | else N=46; |
---|
1161 | } |
---|
1162 | else // Kr (36) |
---|
1163 | { |
---|
1164 | if (rnd<.57) N=48; |
---|
1165 | else if(rnd<.743) N=50; |
---|
1166 | else if(rnd<.859) N=46; |
---|
1167 | else if(rnd<.974) N=47; |
---|
1168 | else if(rnd<.9965) N=44; |
---|
1169 | else N=42; |
---|
1170 | } |
---|
1171 | } |
---|
1172 | } |
---|
1173 | else // ...... Rb - Mo |
---|
1174 | { |
---|
1175 | if (i<40) |
---|
1176 | { |
---|
1177 | if(i==37) // Rb |
---|
1178 | { |
---|
1179 | if (rnd<.7217) N=48; |
---|
1180 | else N=50; |
---|
1181 | } |
---|
1182 | else // SR (38) |
---|
1183 | { |
---|
1184 | if (rnd<.8258) N=50; |
---|
1185 | else if(rnd<.9244) N=48; |
---|
1186 | else if(rnd<.9944) N=49; |
---|
1187 | else N=46; |
---|
1188 | } |
---|
1189 | } |
---|
1190 | else |
---|
1191 | { |
---|
1192 | if(i==40) // Zr |
---|
1193 | { |
---|
1194 | if (rnd<.5145) N=50; |
---|
1195 | else if(rnd<.6883) N=54; |
---|
1196 | else if(rnd<.8598) N=53; |
---|
1197 | else if(rnd<.972) N=51; |
---|
1198 | else N=56; |
---|
1199 | } |
---|
1200 | else // Mo (42) |
---|
1201 | { |
---|
1202 | if (rnd<.2413) N=56; |
---|
1203 | else if(rnd<.4081) N=54; |
---|
1204 | else if(rnd<.5673) N=53; |
---|
1205 | else if(rnd<.7157) N=50; |
---|
1206 | else if(rnd<.8120) N=58; |
---|
1207 | else if(rnd<.9075) N=55; |
---|
1208 | else N=52; |
---|
1209 | } |
---|
1210 | } |
---|
1211 | } |
---|
1212 | } |
---|
1213 | } |
---|
1214 | } |
---|
1215 | else // ====== Ru - U |
---|
1216 | { |
---|
1217 | if (i<66) // ------ Ru - Gd |
---|
1218 | { |
---|
1219 | if (i<54) // ______ Ru - Te |
---|
1220 | { |
---|
1221 | if (i<49) // ...... Ru - Cd |
---|
1222 | { |
---|
1223 | if (i<47) |
---|
1224 | { |
---|
1225 | if(i==44) // Ru |
---|
1226 | { |
---|
1227 | if (rnd<.316) N=58; |
---|
1228 | else if(rnd<.502) N=60; |
---|
1229 | else if(rnd<.673) N=57; |
---|
1230 | else if(rnd<.8) N=55; |
---|
1231 | else if(rnd<.926) N=56; |
---|
1232 | else if(rnd<.9814) N=52; |
---|
1233 | else N=54; |
---|
1234 | } |
---|
1235 | else // Pd (46) |
---|
1236 | { |
---|
1237 | if (rnd<.2733) N=60; |
---|
1238 | else if(rnd<.5379) N=62; |
---|
1239 | else if(rnd<.7612) N=59; |
---|
1240 | else if(rnd<.8784) N=55; |
---|
1241 | else if(rnd<.9898) N=58; |
---|
1242 | else N=56; |
---|
1243 | } |
---|
1244 | } |
---|
1245 | else |
---|
1246 | { |
---|
1247 | if(i==47) // Ag |
---|
1248 | { |
---|
1249 | if(rnd<.51839) N=60; |
---|
1250 | else N=62; |
---|
1251 | } |
---|
1252 | else // Cd (48) |
---|
1253 | { |
---|
1254 | if (rnd<.2873) N=66; |
---|
1255 | else if(rnd<.5286) N=64; |
---|
1256 | else if(rnd<.6566) N=59; |
---|
1257 | else if(rnd<.7815) N=62; |
---|
1258 | else if(rnd<.9037) N=65; |
---|
1259 | else if(rnd<.9786) N=68; |
---|
1260 | else if(rnd<.9911) N=58; |
---|
1261 | else N=60; |
---|
1262 | } |
---|
1263 | } |
---|
1264 | } |
---|
1265 | else // ...... In - Te |
---|
1266 | { |
---|
1267 | if (i<51) |
---|
1268 | { |
---|
1269 | if(i==49) // In |
---|
1270 | { |
---|
1271 | if (rnd<.9577) N=66; |
---|
1272 | else N=64; |
---|
1273 | } |
---|
1274 | else // Sn (50) |
---|
1275 | { |
---|
1276 | if (rnd<.3259) N=70; |
---|
1277 | else if(rnd<.5681) N=68; |
---|
1278 | else if(rnd<.7134) N=66; |
---|
1279 | else if(rnd<.7992) N=69; |
---|
1280 | else if(rnd<.8760) N=67; |
---|
1281 | else if(rnd<.9339) N=74; |
---|
1282 | else if(rnd<.9802) N=72; |
---|
1283 | else if(rnd<.9899) N=62; |
---|
1284 | else N=64; |
---|
1285 | //else if(rnd<.9964) N=64; |
---|
1286 | //else N=65; |
---|
1287 | } |
---|
1288 | } |
---|
1289 | else |
---|
1290 | { |
---|
1291 | if(i==51) // Sb |
---|
1292 | { |
---|
1293 | if (rnd<.5736) N=70; |
---|
1294 | else N=72; |
---|
1295 | } |
---|
1296 | else // Te (52) |
---|
1297 | { |
---|
1298 | if (rnd<.3387) N=78; |
---|
1299 | else if(rnd<.6557) N=76; |
---|
1300 | else if(rnd<.8450) N=74; |
---|
1301 | else if(rnd<.9162) N=73; |
---|
1302 | else if(rnd<.9641) N=72; |
---|
1303 | else if(rnd<.9900) N=70; |
---|
1304 | else if(rnd<.99905) N=71; |
---|
1305 | else N=68; |
---|
1306 | } |
---|
1307 | } |
---|
1308 | } |
---|
1309 | } |
---|
1310 | else // ______ Xe - Gd |
---|
1311 | { |
---|
1312 | if (i<60) // ...... Xe - B |
---|
1313 | { |
---|
1314 | if (i<57) |
---|
1315 | { |
---|
1316 | if(i==54) // Xe |
---|
1317 | { |
---|
1318 | if (rnd<.269) N=78; |
---|
1319 | else if(rnd<.533) N=75; |
---|
1320 | else if(rnd<.745) N=77; |
---|
1321 | else if(rnd<.849) N=80; |
---|
1322 | else if(rnd<.938) N=82; |
---|
1323 | else if(rnd<.979) N=76; |
---|
1324 | else if(rnd<.9981) N=74; |
---|
1325 | else if(rnd<.9991) N=70; |
---|
1326 | else N=72; |
---|
1327 | } |
---|
1328 | else // Ba (56) |
---|
1329 | { |
---|
1330 | if (rnd<.717) N=82; |
---|
1331 | else if(rnd<.8293) N=81; |
---|
1332 | else if(rnd<.9078) N=80; |
---|
1333 | else if(rnd<.97373) N=79; |
---|
1334 | else if(rnd<.99793) N=78; |
---|
1335 | else if(rnd<.99899) N=74; |
---|
1336 | else N=76; |
---|
1337 | } |
---|
1338 | } |
---|
1339 | else |
---|
1340 | { |
---|
1341 | if(i==57) // La |
---|
1342 | { |
---|
1343 | if (rnd<.999098)N=82; |
---|
1344 | else N=81; |
---|
1345 | } |
---|
1346 | else // Ce (58) |
---|
1347 | { |
---|
1348 | if (rnd<.8843) N=82; |
---|
1349 | else if(rnd<.9956) N=84; |
---|
1350 | else if(rnd<.9981) N=80; |
---|
1351 | else N=78; |
---|
1352 | } |
---|
1353 | } |
---|
1354 | } |
---|
1355 | else // ...... Nd - Gd |
---|
1356 | { |
---|
1357 | if (i<63) |
---|
1358 | { |
---|
1359 | if(i==60) // Nd |
---|
1360 | { |
---|
1361 | if (rnd<.2713) N=82; |
---|
1362 | else if(rnd<.5093) N=84; |
---|
1363 | else if(rnd<.6812) N=86; |
---|
1364 | else if(rnd<.8030) N=83; |
---|
1365 | else if(rnd<.8860) N=85; |
---|
1366 | else if(rnd<.9436) N=88; |
---|
1367 | else N=90; |
---|
1368 | } |
---|
1369 | else // Sm (62) |
---|
1370 | { |
---|
1371 | if (rnd<.267) N=90; |
---|
1372 | else if(rnd<.494) N=92; |
---|
1373 | else if(rnd<.644) N=85; |
---|
1374 | else if(rnd<.782) N=87; |
---|
1375 | else if(rnd<.895) N=86; |
---|
1376 | else if(rnd<.969) N=88; |
---|
1377 | else N=82; |
---|
1378 | } |
---|
1379 | } |
---|
1380 | else |
---|
1381 | { |
---|
1382 | if(i==63) // Eu |
---|
1383 | { |
---|
1384 | if (rnd<.522) N=90; |
---|
1385 | else N=89; |
---|
1386 | } |
---|
1387 | else // Gd (64) |
---|
1388 | { |
---|
1389 | if (rnd<.2484) N=94; |
---|
1390 | else if(rnd<.4670) N=96; |
---|
1391 | else if(rnd<.6717) N=92; |
---|
1392 | else if(rnd<.8282) N=93; |
---|
1393 | else if(rnd<.9762) N=91; |
---|
1394 | else if(rnd<.9980) N=90; |
---|
1395 | else N=88; |
---|
1396 | } |
---|
1397 | } |
---|
1398 | } |
---|
1399 | } |
---|
1400 | } |
---|
1401 | else // ------ Dy - U |
---|
1402 | { |
---|
1403 | if (i<76) // ______ Dy - Re |
---|
1404 | { |
---|
1405 | if (i<72) // ...... Dy - Lu |
---|
1406 | { |
---|
1407 | if (i<70) |
---|
1408 | { |
---|
1409 | if(i==66) // Dy |
---|
1410 | { |
---|
1411 | if (rnd<.282) N=98; |
---|
1412 | else if(rnd<.537) N=96; |
---|
1413 | else if(rnd<.786) N=97; |
---|
1414 | else if(rnd<.975) N=95; |
---|
1415 | else if(rnd<.9984) N=94; |
---|
1416 | else if(rnd<.9994) N=92; |
---|
1417 | else N=90; |
---|
1418 | } |
---|
1419 | else // Er (68) |
---|
1420 | { |
---|
1421 | if (rnd<.3360) N= 98; |
---|
1422 | else if(rnd<.6040) N=100; |
---|
1423 | else if(rnd<.8335) N= 99; |
---|
1424 | else if(rnd<.9825) N=102; |
---|
1425 | else if(rnd<.9986) N= 96; |
---|
1426 | else N= 94; |
---|
1427 | } |
---|
1428 | } |
---|
1429 | else |
---|
1430 | { |
---|
1431 | if(i==70) // Yb |
---|
1432 | { |
---|
1433 | if (rnd<.3180) N=104; |
---|
1434 | else if(rnd<.5370) N=102; |
---|
1435 | else if(rnd<.6982) N=103; |
---|
1436 | else if(rnd<.8412) N=101; |
---|
1437 | else if(rnd<.9682) N=106; |
---|
1438 | else if(rnd<.9987) N=100; |
---|
1439 | else N= 98; |
---|
1440 | } |
---|
1441 | else // Lu (71) |
---|
1442 | { |
---|
1443 | if (rnd<.9741) N=104; |
---|
1444 | else N=105; |
---|
1445 | } |
---|
1446 | } |
---|
1447 | } |
---|
1448 | else // ...... Hf - Re |
---|
1449 | { |
---|
1450 | if (i<74) |
---|
1451 | { |
---|
1452 | if(i==72) // Hf |
---|
1453 | { |
---|
1454 | if (rnd<.35100) N=108; |
---|
1455 | else if(rnd<.62397) N=106; |
---|
1456 | else if(rnd<.81003) N=105; |
---|
1457 | else if(rnd<.94632) N=107; |
---|
1458 | else if(rnd<.99838) N=104; |
---|
1459 | else N=102; |
---|
1460 | } |
---|
1461 | else // Ta (73) |
---|
1462 | { |
---|
1463 | if(rnd<.99988) N=108; |
---|
1464 | else N=107; |
---|
1465 | } |
---|
1466 | } |
---|
1467 | else |
---|
1468 | { |
---|
1469 | if(i==74) // W |
---|
1470 | { |
---|
1471 | if (rnd<.307) N=110; |
---|
1472 | else if(rnd<.593) N=112; |
---|
1473 | else if(rnd<.856) N=108; |
---|
1474 | else if(rnd<.9988) N=109; |
---|
1475 | else N=106; |
---|
1476 | } |
---|
1477 | else // Re (75) |
---|
1478 | { |
---|
1479 | if (rnd<.626) N=112; |
---|
1480 | else N=110; |
---|
1481 | } |
---|
1482 | } |
---|
1483 | } |
---|
1484 | } |
---|
1485 | else // ______ Os - U |
---|
1486 | { |
---|
1487 | if (i<81) // ...... Os - Hg |
---|
1488 | { |
---|
1489 | if (i<78) |
---|
1490 | { |
---|
1491 | if(i==76) // Os |
---|
1492 | { |
---|
1493 | if (rnd<.410) N=116; |
---|
1494 | else if(rnd<.674) N=114; |
---|
1495 | else if(rnd<.835) N=113; |
---|
1496 | else if(rnd<.968) N=112; |
---|
1497 | else if(rnd<.984) N=111; |
---|
1498 | else if(rnd<.9998) N=110; |
---|
1499 | else N=108; |
---|
1500 | } |
---|
1501 | else // Ir (77) |
---|
1502 | { |
---|
1503 | if (rnd<.627) N=116; |
---|
1504 | else N=114; |
---|
1505 | } |
---|
1506 | } |
---|
1507 | else |
---|
1508 | { |
---|
1509 | if(i==78) // Pt |
---|
1510 | { |
---|
1511 | if (rnd<.338) N=117; |
---|
1512 | else if(rnd<.667) N=116; |
---|
1513 | else if(rnd<.920) N=118; |
---|
1514 | else if(rnd<.992) N=120; |
---|
1515 | else if(rnd<.9999) N=114; |
---|
1516 | else N=112; |
---|
1517 | } |
---|
1518 | else // Hg (80) |
---|
1519 | { |
---|
1520 | if (rnd<.2986) N=122; |
---|
1521 | else if(rnd<.5296) N=120; |
---|
1522 | else if(rnd<.6983) N=119; |
---|
1523 | else if(rnd<.8301) N=121; |
---|
1524 | else if(rnd<.9298) N=118; |
---|
1525 | else if(rnd<.9985) N=124; |
---|
1526 | else N=116; |
---|
1527 | } |
---|
1528 | } |
---|
1529 | } |
---|
1530 | else // ...... Tl - U |
---|
1531 | { |
---|
1532 | if (i<92) |
---|
1533 | { |
---|
1534 | if (i==81) // Tl |
---|
1535 | { |
---|
1536 | if (rnd<.70476) N=124; |
---|
1537 | else N=122; |
---|
1538 | } |
---|
1539 | else // Pb (82) |
---|
1540 | { |
---|
1541 | if (rnd<.524) N=126; |
---|
1542 | else if(rnd<.765) N=124; |
---|
1543 | else if(rnd<.986) N=125; |
---|
1544 | else N=122; |
---|
1545 | } |
---|
1546 | } |
---|
1547 | else // U (92) |
---|
1548 | { |
---|
1549 | if (rnd<.992745)N=146; |
---|
1550 | else if(rnd<.999945)N=143; |
---|
1551 | else N=142; |
---|
1552 | } |
---|
1553 | } |
---|
1554 | } |
---|
1555 | } |
---|
1556 | } |
---|
1557 | } |
---|
1558 | else if(N<0) |
---|
1559 | { |
---|
1560 | N=-N; |
---|
1561 | G4cerr<<"---Worn---G4QIsotope::RandomizeNeutrons:UnstableElement's used Z="<<i<<G4endl; |
---|
1562 | } |
---|
1563 | #ifdef debug |
---|
1564 | G4cout<<"G4QIsotope::RandomizeNeutrons: Z="<<i<<", N="<<N<<G4endl; |
---|
1565 | #endif |
---|
1566 | return N; |
---|
1567 | } |
---|
1568 | |
---|
1569 | // Returns the input index (if it is >0 & unique) or theFirstFreeIndex (<=0 or nonunique) |
---|
1570 | G4int G4QIsotope::InitElement(G4int Z, G4int index, // Ret: -1 - Empty, -2 - Wrong (sum>1) |
---|
1571 | vector<pair<G4int,G4double>*>* abund) |
---|
1572 | // ======================================================================= |
---|
1573 | { |
---|
1574 | G4int I=abund->size(); |
---|
1575 | #ifdef debug |
---|
1576 | G4cout<<"G4QIsotope::InitElement: called with I="<<I<<" pairs,Z="<<Z<<",i="<<indexG4endl; |
---|
1577 | #endif |
---|
1578 | if(I<=0) |
---|
1579 | { |
---|
1580 | G4cerr<<"--Worning--G4QIsotope::InitEl:(-1)0VectorOfNewEl,Z="<<Z<<",i="<<index<<G4endl; |
---|
1581 | return -2; |
---|
1582 | } |
---|
1583 | if(IsDefined(Z,index)) // This index is already defined |
---|
1584 | { |
---|
1585 | G4cerr<<"-Worning-G4QIsotope::InitEl:VONewEl,Z="<<Z<<",ind="<<index<<" exists"<<G4endl; |
---|
1586 | return index; |
---|
1587 | } |
---|
1588 | G4int ZInd=1000*index+Z; // Fake Z increased by the UserDefinedIndex |
---|
1589 | vector<pair<G4int,G4double>*>*A = new vector<pair<G4int,G4double>*>; |
---|
1590 | vector<pair<G4int,G4double>*>*S = new vector<pair<G4int,G4double>*>; |
---|
1591 | vector<pair<G4int,G4double>*>*C = new vector<pair<G4int,G4double>*>; |
---|
1592 | #ifdef debug |
---|
1593 | G4cout<<"G4QIsotope::InitElement: A & S & C vectors are alocated"<<G4endl; |
---|
1594 | #endif |
---|
1595 | G4double sumAbu=0; // Summ of abbundancies |
---|
1596 | for(G4int j=0; j<I; j++) |
---|
1597 | { |
---|
1598 | G4int N=(*abund)[j]->first; |
---|
1599 | G4double abu=(*abund)[j]->second; |
---|
1600 | #ifdef debug |
---|
1601 | G4cout<<"G4QIsotope::InitElement: pair#"<<j<<", N="<<N<<", abund="<<abu<<G4endl; |
---|
1602 | #endif |
---|
1603 | sumAbu+=abu; |
---|
1604 | if(j==I-1.) |
---|
1605 | { |
---|
1606 | if(fabs(sumAbu-1.)>.00001) |
---|
1607 | { |
---|
1608 | G4cerr<<"--Worning--G4QIsotope::InitEl:maxSum="<<sumAbu<<" is fixed to 1."<<G4endl; |
---|
1609 | abu+=1.-sumAbu; |
---|
1610 | sumAbu=1.; |
---|
1611 | } |
---|
1612 | else if(sumAbu-abu>1.) |
---|
1613 | { |
---|
1614 | G4cerr<<"--Worning--G4QIsotope::InitEl:(-2)WrongAbund,Z="<<Z<<",i="<<index<<G4endl; |
---|
1615 | for(G4int k=0; k<I-1; k++) |
---|
1616 | { |
---|
1617 | delete (*A)[k]; |
---|
1618 | delete (*S)[k]; |
---|
1619 | delete (*C)[k]; |
---|
1620 | } |
---|
1621 | delete A; |
---|
1622 | delete S; |
---|
1623 | delete C; |
---|
1624 | return -2; |
---|
1625 | } |
---|
1626 | #ifdef debug |
---|
1627 | G4cout<<"G4QIsotope::InitElement:TheLastIsChecked it isOK or coredTo "<<sumAbu<<G4endl; |
---|
1628 | #endif |
---|
1629 | } |
---|
1630 | pair<G4int,G4double>* abP= new pair<G4int,G4double>(N,abu); |
---|
1631 | A->push_back(abP); |
---|
1632 | pair<G4int,G4double>* saP= new pair<G4int,G4double>(N,sumAbu); |
---|
1633 | S->push_back(saP); |
---|
1634 | pair<G4int,G4double>* csP= new pair<G4int,G4double>(N,0.); |
---|
1635 | C->push_back(csP); |
---|
1636 | #ifdef debug |
---|
1637 | G4cout<<"G4QIsotope::InitElement: A & S & C are filled nP="<<C->size()<<G4endl; |
---|
1638 | #endif |
---|
1639 | } |
---|
1640 | pair<G4int,vector<pair<G4int,G4double>*>*>* newAP= |
---|
1641 | new pair<G4int,vector<pair<G4int,G4double>*>*>(ZInd,A); |
---|
1642 | newElems.push_back(newAP); |
---|
1643 | pair<G4int,vector<pair<G4int,G4double>*>*>* newSA= |
---|
1644 | new pair<G4int,vector<pair<G4int,G4double>*>*>(ZInd,S); |
---|
1645 | newSumAb.push_back(newSA); |
---|
1646 | pair<G4int,vector<pair<G4int,G4double>*>*>* newCP= |
---|
1647 | new pair<G4int,vector<pair<G4int,G4double>*>*>(ZInd,C); |
---|
1648 | newIsoCS.push_back(newCP); |
---|
1649 | #ifdef debug |
---|
1650 | G4cout<<"G4QIsotope::InitElement: newElems & newSumAb & newIsoCS are filled "<<G4endl; |
---|
1651 | #endif |
---|
1652 | return index; |
---|
1653 | } |
---|
1654 | |
---|
1655 | // The highest index defined for Element with Z (Index>0 correspondToUserDefinedElements) |
---|
1656 | G4int G4QIsotope::GetLastIndex(G4int Z) // Returns theLastDefinedIndex (if onlyNatural: =0) |
---|
1657 | // ================================= |
---|
1658 | { |
---|
1659 | #ifdef debug |
---|
1660 | G4cout<<"G4QIsotope::GetLastIndex is called Z="<<Z<<G4endl; |
---|
1661 | #endif |
---|
1662 | G4int mind=0; // Prototype of the maximum existing index for this Z |
---|
1663 | G4int nE=newElems.size(); // A number of definitions in the newElements vector |
---|
1664 | if(nE) for(G4int i=0; i<nE; i++) // LOOP over new UserDefinedElements |
---|
1665 | { |
---|
1666 | G4int zin=newElems[i]->first; |
---|
1667 | G4int zi=zin%1000; // Existing Z |
---|
1668 | G4int in=zin/1000; // Existing index |
---|
1669 | if(Z==zi && in>mind) mind=in; // maximum index for this Z |
---|
1670 | } |
---|
1671 | return mind; |
---|
1672 | } |
---|
1673 | |
---|
1674 | // Indices can have differen numbers (not 1,2,3,...) & in different sequences (9,3,7,...) |
---|
1675 | G4bool G4QIsotope::IsDefined(G4int Z, G4int Ind) // Ind is an index to be found (true) |
---|
1676 | // ========================================== |
---|
1677 | { |
---|
1678 | #ifdef debug |
---|
1679 | G4cout<<"G4QIsotope::IsDefined is called Z="<<Z<<", I="<<Ind<<G4endl; |
---|
1680 | #endif |
---|
1681 | if(Ind<=0) |
---|
1682 | { |
---|
1683 | if(Ind<0) G4cerr<<"-W-G4QIsotope::IsDefined: Z="<<Z<<", Ind="<<Ind<<" < 0->=0"<<G4endl; |
---|
1684 | return true; // to avoid definition with the negative index |
---|
1685 | } |
---|
1686 | G4int nE=newElems.size(); // A number of definitions in the newElements vector |
---|
1687 | if(nE) for(G4int i=0; i<nE; i++) // LOOP over new UserDefinedElements |
---|
1688 | { |
---|
1689 | G4int zin=newElems[i]->first; |
---|
1690 | G4int zi=zin%1000; // Existing Z |
---|
1691 | G4int in=zin/1000; // Existing index |
---|
1692 | if(Z==zi && Ind==in) return true; // The index for the element Z is found |
---|
1693 | } |
---|
1694 | return false; // The index for the element Z is not found |
---|
1695 | } |
---|
1696 | |
---|
1697 | // A#ofNeutrons in theElement with Z & UserDefIndex. Universal for Nat(index=0) & UserDefEl |
---|
1698 | G4int G4QIsotope::GetNeutrons(G4int Z, G4int index) // If theElem doesn't exist, returns <0 |
---|
1699 | // ============================================= |
---|
1700 | { |
---|
1701 | #ifdef debug |
---|
1702 | G4cout<<"G4QIsotope::GetNeutrons is called Z="<<Z<<", index="<<index<<G4endl; |
---|
1703 | #endif |
---|
1704 | // To reduce the code, but make the member function a bit slower, one can use for natural |
---|
1705 | // isotopes the same algorithm as for the newElements, splitting the natElements Vector |
---|
1706 | if(!index) return RandomizeNeutrons(Z); // @@ Fast decision for the natural isotopes |
---|
1707 | else if(index<0) |
---|
1708 | { |
---|
1709 | G4cerr<<"---Worning---G4QIsotope::GetNeutrons:(-2) Negative Index i="<<index<<G4endl; |
---|
1710 | return -2; |
---|
1711 | } |
---|
1712 | // For the positive index tries to randomize the newUserDefinedElement |
---|
1713 | G4bool found=false; // Prototype of the"ZWithTheSameIndex is found" event |
---|
1714 | G4int nE=newElems.size(); // A number of definitions in the newElements Vector |
---|
1715 | G4int i=0; |
---|
1716 | if(nE) for(i=0; i<nE; i++) |
---|
1717 | { |
---|
1718 | G4int zin=newElems[i]->first; |
---|
1719 | G4int in=zin/1000; // Existing index |
---|
1720 | G4int zi=zin%1000; // Existing Z |
---|
1721 | if(Z==zi && in==index) |
---|
1722 | { |
---|
1723 | found=true; // The newElement with the same Z & index is found |
---|
1724 | break; // Finish the search and quit the loop |
---|
1725 | } |
---|
1726 | } |
---|
1727 | if(!found) |
---|
1728 | { |
---|
1729 | G4cerr<<"--Worning--G4QIsotope::GetNeutrons:(-1) NotFound Z="<<Z<<",i="<<index<<G4endl; |
---|
1730 | return -1; |
---|
1731 | } |
---|
1732 | vector<pair<G4int,G4double>*>* abu = newSumAb[i]->second; |
---|
1733 | G4int nn = abu->size(); // A#Of UserDefinedIsotopes for the newElement |
---|
1734 | if(nn>0) |
---|
1735 | { |
---|
1736 | if(nn==1) return (*abu)[0]->first; |
---|
1737 | else |
---|
1738 | { |
---|
1739 | G4double rnd=G4UniformRand(); |
---|
1740 | G4int j=0; |
---|
1741 | for(j=0; j<nn; j++) if ( rnd < (*abu)[j]->second ) break; |
---|
1742 | if(j>=nn) j=nn-1; |
---|
1743 | return (*abu)[j]->first; |
---|
1744 | } |
---|
1745 | } |
---|
1746 | else |
---|
1747 | { |
---|
1748 | G4cerr<<"--Worning--G4QIsotope::GetNeutrons:(-3) Empty Z="<<Z<<",i="<<index<<G4endl; |
---|
1749 | return -3; |
---|
1750 | } |
---|
1751 | } |
---|
1752 | |
---|
1753 | // Get a pointer to the vector of pairs(N,CrosSec), where N is used to calculate CrosSec |
---|
1754 | vector<pair<G4int,G4double>*>* G4QIsotope::GetCSVector(G4int Z, G4int index) |
---|
1755 | // ============================================= |
---|
1756 | { |
---|
1757 | #ifdef debug |
---|
1758 | G4cout<<"G4QIsotope::GetCSVector is called"<<G4endl; |
---|
1759 | #endif |
---|
1760 | if(index<0) |
---|
1761 | { |
---|
1762 | G4cerr<<"---Worning---G4QIsotope::GetSCVector:(-1) Negative Index i="<<index<<G4endl; |
---|
1763 | return 0; |
---|
1764 | } |
---|
1765 | else if(!index) return natIsoCrosS[Z]; |
---|
1766 | // For the positive index tries to find the newUserDefinedElement |
---|
1767 | G4bool found=false; // Prototype of the"ZWithTheSameIndex is found" event |
---|
1768 | G4int nE=newIsoCS.size(); // A number of definitions in the newElements Vector |
---|
1769 | G4int i=0; |
---|
1770 | if(nE) for(i=0; i<nE; i++) |
---|
1771 | { |
---|
1772 | G4int zin=newIsoCS[i]->first; |
---|
1773 | G4int in=zin/1000; // Existing index |
---|
1774 | G4int zi=zin%1000; // Existing Z |
---|
1775 | if(Z==zi && in==index) |
---|
1776 | { |
---|
1777 | found=true; // The newElement with the same Z & index is found |
---|
1778 | break; // Finish the search and quit the loop |
---|
1779 | } |
---|
1780 | } |
---|
1781 | if(!found) |
---|
1782 | { |
---|
1783 | G4cerr<<"--Worning--G4QIsotope::GetSCVector:(-2) NotFound Z="<<Z<<",i="<<index<<G4endl; |
---|
1784 | return 0; |
---|
1785 | } |
---|
1786 | return newIsoCS[i]->second; |
---|
1787 | } |
---|
1788 | |
---|
1789 | // Get a pointer to the vector of pairs(N,IntAbundancy) for the element with Z |
---|
1790 | vector<pair<G4int,G4double>*>* G4QIsotope::GetAbuVector(G4int Z, G4int index) |
---|
1791 | // ============================================== |
---|
1792 | { |
---|
1793 | #ifdef debug |
---|
1794 | G4cout<<"G4QIsotope::GetAbuVector is called"<<G4endl; |
---|
1795 | #endif |
---|
1796 | if(index<0) |
---|
1797 | { |
---|
1798 | G4cerr<<"---Worning---G4QIsotope::GetAbuVector:(-1) Negative Index i="<<index<<G4endl; |
---|
1799 | return 0; |
---|
1800 | } |
---|
1801 | else if(!index) return natElements[Z]; |
---|
1802 | // For the positive index tries to find the newUserDefinedElement |
---|
1803 | G4bool found=false; // Prototype of the"ZWithTheSameIndex is found" event |
---|
1804 | G4int nE=newElems.size(); // A number of definitions in the newElements Vector |
---|
1805 | G4int i=0; |
---|
1806 | if(nE) for(i=0; i<nE; i++) |
---|
1807 | { |
---|
1808 | G4int zin=newElems[i]->first; |
---|
1809 | G4int in=zin/1000; // Existing index |
---|
1810 | G4int zi=zin%1000; // Existing Z |
---|
1811 | if(Z==zi && in==index) |
---|
1812 | { |
---|
1813 | found=true; // The newElement with the same Z & index is found |
---|
1814 | break; // Finish the search and quit the loop |
---|
1815 | } |
---|
1816 | } |
---|
1817 | if(!found) |
---|
1818 | { |
---|
1819 | G4cerr<<"--Worning--G4QIsotope::GetAbuVector:(-2)NotFound Z="<<Z<<",i="<<index<<G4endl; |
---|
1820 | return 0; |
---|
1821 | } |
---|
1822 | return newElems[i]->second; |
---|
1823 | } |
---|
1824 | |
---|
1825 | // Get a pointer to the vector of pairs(N,SumAbundancy) for the element with Z |
---|
1826 | vector<pair<G4int,G4double>*>* G4QIsotope::GetSumAVector(G4int Z, G4int index) |
---|
1827 | // =============================================== |
---|
1828 | { |
---|
1829 | #ifdef debug |
---|
1830 | G4cout<<"G4QIsotope::GetSumAVector is called"<<G4endl; |
---|
1831 | #endif |
---|
1832 | if(index<0) |
---|
1833 | { |
---|
1834 | G4cerr<<"---Worning---G4QIsotope::GetSumAVector:(-1) Negative Index i="<<index<<G4endl; |
---|
1835 | return 0; |
---|
1836 | } |
---|
1837 | else if(!index) return natSumAbund[Z]; |
---|
1838 | // For the positive index tries to find the newUserDefinedElement |
---|
1839 | G4bool found=false; // Prototype of the"ZWithTheSameIndex is found" event |
---|
1840 | G4int nE=newSumAb.size(); // A number of definitions in the newElements Vector |
---|
1841 | G4int i=0; |
---|
1842 | if(nE) for(i=0; i<nE; i++) |
---|
1843 | { |
---|
1844 | G4int zin=newSumAb[i]->first; |
---|
1845 | G4int in=zin/1000; // Existing index |
---|
1846 | G4int zi=zin%1000; // Existing Z |
---|
1847 | if(Z==zi && in==index) |
---|
1848 | { |
---|
1849 | found=true; // The newElement with the same Z & index is found |
---|
1850 | break; // Finish the search and quit the loop |
---|
1851 | } |
---|
1852 | } |
---|
1853 | if(!found) |
---|
1854 | { |
---|
1855 | G4cerr<<"-Worning-G4QIsotope::GetSumAVector:(-2)Not Found Z="<<Z<<",i="<<index<<G4endl; |
---|
1856 | return 0; |
---|
1857 | } |
---|
1858 | return newSumAb[i]->second; |
---|
1859 | } |
---|
1860 | |
---|
1861 | // Calculates the mean Cross Section for the initialized Element(ind=0 Nat,ind>0 UserDef) |
---|
1862 | G4double G4QIsotope::GetMeanCrossSection(G4int Z, G4int index) |
---|
1863 | // ========================================= |
---|
1864 | { |
---|
1865 | vector<pair<G4int,G4double>*>* ab; |
---|
1866 | vector<pair<G4int,G4double>*>* cs; |
---|
1867 | #ifdef ppdebug |
---|
1868 | G4cout<<"G4QIsotope::GetMeanCrossSection is called"<<G4endl; |
---|
1869 | #endif |
---|
1870 | if(index<0) |
---|
1871 | { |
---|
1872 | G4cerr<<"---Worning---G4QIsotope::GetMeanCS:(-1) Negative Index i="<<index<<G4endl; |
---|
1873 | return -1.; |
---|
1874 | } |
---|
1875 | else if(!index) // =========> Natural Abundancies for Isotopes of the Element |
---|
1876 | { |
---|
1877 | #ifdef ppdebug |
---|
1878 | G4cout<<"G4QIsotope::GetMeanCrossSection: Nat Abundance, Z="<<Z<<G4endl; |
---|
1879 | #endif |
---|
1880 | ab=natElements[Z]; |
---|
1881 | cs=natIsoCrosS[Z]; |
---|
1882 | } |
---|
1883 | else // =========> UserDefinedAbundancies for Isotopes of theElement |
---|
1884 | { |
---|
1885 | #ifdef ppdebug |
---|
1886 | G4cout<<"G4QIsotope::GetMeanCrossSection: Art Abund, Z="<<Z<<",ind="<<index<<G4endl; |
---|
1887 | #endif |
---|
1888 | // For the positive index tries to find the newUserDefinedElement |
---|
1889 | G4bool found=false; // Prototype of the"ZWithTheSameIndex is found" event |
---|
1890 | G4int nE=newIsoCS.size(); // A number of definitions in the newElements Vector |
---|
1891 | G4int i=0; |
---|
1892 | if(nE) for(i=0; i<nE; i++) |
---|
1893 | { |
---|
1894 | G4int zin=newIsoCS[i]->first; |
---|
1895 | G4int in=zin/1000; // Existing index |
---|
1896 | G4int zi=zin%1000; // Existing Z |
---|
1897 | if(Z==zi && in==index) |
---|
1898 | { |
---|
1899 | found=true; // The newElement with the same Z & index is found |
---|
1900 | break; // Finish the search and quit the loop |
---|
1901 | } |
---|
1902 | } |
---|
1903 | if(!found) |
---|
1904 | { |
---|
1905 | G4cerr<<"--Worning--G4QIsotope::GetMeanCS:(-2) NotFound Z="<<Z<<",i="<<index<<G4endl; |
---|
1906 | return -2.; |
---|
1907 | } |
---|
1908 | ab=newElems[i]->second; |
---|
1909 | cs=newIsoCS[i]->second; |
---|
1910 | } |
---|
1911 | G4int nis=ab->size(); |
---|
1912 | //G4double last=0.; |
---|
1913 | if(!nis) |
---|
1914 | { |
---|
1915 | G4cerr<<"--Worning--G4QIsotope::GetMeanCS:(-3) Empty Z="<<Z<<",i="<<index<<G4endl; |
---|
1916 | return -3.; |
---|
1917 | } |
---|
1918 | else |
---|
1919 | { |
---|
1920 | G4double sum=0.; |
---|
1921 | for(G4int j=0; j<nis; j++) |
---|
1922 | { |
---|
1923 | G4double cur=(*ab)[j]->second; |
---|
1924 | //G4double abunda=cur-last; |
---|
1925 | //last=cur; |
---|
1926 | #ifdef ppdebug |
---|
1927 | G4cout<<"G4QIsot::GetMeanCS:j="<<j<<",ab="<<cur<<",CS="<<(*cs)[j]->second<<G4endl; |
---|
1928 | #endif |
---|
1929 | //sum+=abunda * (*cs)[j]->second; |
---|
1930 | sum+=cur * (*cs)[j]->second; |
---|
1931 | } |
---|
1932 | return sum; |
---|
1933 | } |
---|
1934 | } |
---|
1935 | |
---|
1936 | // Randomize A#OfNeutrons in the Isotope weighted by theAbubdancies and theCrossSections |
---|
1937 | G4int G4QIsotope::GetCSNeutrons(G4int Z, G4int index) |
---|
1938 | // =================================== |
---|
1939 | { |
---|
1940 | vector<pair<G4int,G4double>*>* ab; |
---|
1941 | vector<pair<G4int,G4double>*>* cs; |
---|
1942 | #ifdef debug |
---|
1943 | G4cout<<"G4QIsotope::GetCSNeutrons is called"<<G4endl; |
---|
1944 | #endif |
---|
1945 | if(index<0) |
---|
1946 | { |
---|
1947 | G4cerr<<"---Worning---G4QIsotope::GetCSNeutrons:(-1) Negative Index i="<<index<<G4endl; |
---|
1948 | return -1; |
---|
1949 | } |
---|
1950 | else if(!index) // =========> Natural Abundancies for Isotopes of the Element |
---|
1951 | { |
---|
1952 | ab=natElements[Z]; |
---|
1953 | cs=natIsoCrosS[Z]; |
---|
1954 | } |
---|
1955 | else // =========> UserDefinedAbundancies for Isotopes of theElement |
---|
1956 | { |
---|
1957 | // For the positive index tries to find the newUserDefinedElement |
---|
1958 | G4bool found=false; // Prototype of the"ZWithTheSameIndex is found" event |
---|
1959 | G4int nE=newIsoCS.size(); // A number of definitions in the newElements Vector |
---|
1960 | G4int i=0; |
---|
1961 | if(nE) for(i=0; i<nE; i++) |
---|
1962 | { |
---|
1963 | G4int zin=newIsoCS[i]->first; |
---|
1964 | G4int in=zin/1000; // Existing index |
---|
1965 | G4int zi=zin%1000; // Existing Z |
---|
1966 | if(Z==zi && in==index) |
---|
1967 | { |
---|
1968 | found=true; // The newElement with the same Z & index is found |
---|
1969 | break; // Finish the search and quit the loop |
---|
1970 | } |
---|
1971 | } |
---|
1972 | if(!found) |
---|
1973 | { |
---|
1974 | G4cerr<<"--Worning--G4QIsotope::GetCSNeut:(-2) NotFound Z="<<Z<<",i="<<index<<G4endl; |
---|
1975 | return -2; |
---|
1976 | } |
---|
1977 | ab=newElems[i]->second; |
---|
1978 | cs=newIsoCS[i]->second; |
---|
1979 | } |
---|
1980 | G4int nis=ab->size(); |
---|
1981 | G4double last=0.; |
---|
1982 | if(!nis) |
---|
1983 | { |
---|
1984 | G4cerr<<"--Worning--G4QIsotope::GetCSNeutrons:(-3) Empty Z="<<Z<<",i="<<index<<G4endl; |
---|
1985 | return -3; |
---|
1986 | } |
---|
1987 | else |
---|
1988 | { |
---|
1989 | G4double sum=0.; |
---|
1990 | vector<G4double> scs(nis); |
---|
1991 | for(G4int j=0; j<nis; j++) |
---|
1992 | { |
---|
1993 | G4double cur=(*ab)[j]->second; |
---|
1994 | G4double abunda=cur-last; |
---|
1995 | last=cur; |
---|
1996 | sum+=abunda * (*cs)[j]->second;; |
---|
1997 | scs.push_back(sum); |
---|
1998 | } |
---|
1999 | G4double rnd=sum*G4UniformRand(); |
---|
2000 | sum=0; |
---|
2001 | G4int k=0; |
---|
2002 | if(nis>1) for(k=0; k<nis; k++) if(rnd<scs[k]) break; |
---|
2003 | return (*ab)[k]->first; |
---|
2004 | } |
---|
2005 | } |
---|