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

Last change on this file since 1340 was 1340, checked in by garnier, 14 years ago

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