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