Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros
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