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