Context Navigation

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

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

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

Download in other formats: