Context Navigation

source: trunk/source/processes/hadronic/models/chiral_inv_phase_space/body/src/G4QIsotope.cc @ 1071

Last change on this file since 1071 was 1055, checked in by garnier, 15 years ago
maj sur la beta de geant 4.9.3
File size: 71.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: