G4QIsotope.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 //      ---------------- G4QIsotope class ----------------
00030 //             by Mikhail Kossov, December 2003.
00031 //  class G4QIsotope gives Isotopes for Elements (CHIPS branch of G4)
00032 // ------------------------------------------------------------------
00033 // ****************************************************************************************
00034 // ********** This CLASS is temporary moved from the photolepton_hadron directory *********
00035 // ******* DO NOT MAKE ANY CHANGE! With time it'll move back to photolepton...(M.K.) ******
00036 // ****************************************************************************************
00037 //
00038 //       1         2         3         4         5         6         7         8         9
00039 //34567890123456789012345678901234567890123456789012345678901234567890123456789012345678901
00040 // ------------------------------------------------------------------------------
00041 // Short description: containes the natural abundance DB for isotops and permits
00042 // new artificial materials with unnatural abundance (e.g. enreached Deuterium).
00043 // Uses isotope cross-sections for calculation of the mean cross-sections for the
00044 // Element (fixed Z).
00045 // ------------------------------------------------------------------------------
00046 
00047 
00048 //#define debug
00049 //#define cdebug
00050 //#define pdebug
00051 //#define ppdebug
00052 //#define sdebug
00053 
00054 #include "G4QIsotope.hh"
00055 #include <cmath>
00056 using namespace std;
00057 
00058 vector<vector<pair<G4int,G4double>*>*>G4QIsotope::natElements;//NaturalElems
00059 vector<vector<pair<G4int,G4double>*>*>G4QIsotope::natSumAbund;//NaturalSumAb
00060 vector<vector<pair<G4int,G4double>*>*>G4QIsotope::natIsoCrosS;//CSOfNatElems
00061 vector<pair<G4int,vector<pair<G4int,G4double>*>*>*>G4QIsotope::newElems;
00062 vector<pair<G4int,vector<pair<G4int,G4double>*>*>*>G4QIsotope::newSumAb;
00063 vector<pair<G4int,vector<pair<G4int,G4double>*>*>*>G4QIsotope::newIsoCS;
00064 
00065 G4QIsotope::G4QIsotope() 
00066 {
00067 #ifdef cdebug
00068   G4cout<<"G4QIsotope::Constructor is called"<<G4endl;
00069 #endif
00070   vector<vector<pair<G4int,G4double> >*> natEl;
00071 #ifdef cdebug
00072   G4cout<<"G4QIsotope::Constructor natEl is booked"<<G4endl;
00073 #endif
00074   vector<pair<G4int,G4double> >*a0=new vector<pair<G4int,G4double> >;
00075 #ifdef cdebug
00076   G4cout<<"G4QIsotope::Constructor a0 is booked"<<G4endl;
00077 #endif
00078   a0->push_back(make_pair(1,1.));
00079 #ifdef cdebug
00080   G4cout<<"G4QIsotope::Constructor a0 is filled by a pair"<<G4endl;
00081 #endif
00082   natEl.push_back(a0);
00083 #ifdef cdebug
00084   G4cout<<"G4QIsotope::Constructor a0 is filled in natEl"<<G4endl;
00085 #endif
00086   // If an error is found in this initialization, please, correct the simple tree below
00087   vector<pair<G4int,G4double> >*a1=new vector<pair<G4int,G4double> >; // 1-H
00088   a1->push_back(make_pair(0,.99985));
00089   a1->push_back(make_pair(1,1.));
00090   natEl.push_back(a1);
00091   vector<pair<G4int,G4double> >*a2=new vector<pair<G4int,G4double> >; // 2-He
00092   a2->push_back(make_pair(2,.999999863));
00093   a2->push_back(make_pair(1,1.));
00094   natEl.push_back(a2);
00095   vector<pair<G4int,G4double> >*a3=new vector<pair<G4int,G4double> >; // 3-Li
00096   a3->push_back(make_pair(4,.925));
00097   a3->push_back(make_pair(3,1.));
00098   natEl.push_back(a3);
00099   vector<pair<G4int,G4double> >*a4=new vector<pair<G4int,G4double> >; // 4-Be
00100   a4->push_back(make_pair(5,1.));
00101   natEl.push_back(a4);
00102   vector<pair<G4int,G4double> >*a5=new vector<pair<G4int,G4double> >; // 5-B
00103   a5->push_back(make_pair(6,.801));
00104   a5->push_back(make_pair(5,1.));
00105   natEl.push_back(a5);
00106   vector<pair<G4int,G4double> >*a6=new vector<pair<G4int,G4double> >; // 6-C
00107   a6->push_back(make_pair(6,.989));
00108   a6->push_back(make_pair(7,1.));
00109   natEl.push_back(a6);
00110   vector<pair<G4int,G4double> >*a7=new vector<pair<G4int,G4double> >; // 7-N
00111   a7->push_back(make_pair(7,.9963));
00112   a7->push_back(make_pair(8,1.));
00113   natEl.push_back(a7);
00114   vector<pair<G4int,G4double> >*a8=new vector<pair<G4int,G4double> >; // 8-O
00115   a8->push_back(make_pair(8,.9976));
00116   a8->push_back(make_pair(10,.9996));
00117   a8->push_back(make_pair(9,1.));
00118   natEl.push_back(a8);
00119   vector<pair<G4int,G4double> >*a9=new vector<pair<G4int,G4double> >; // 9-F
00120   a9->push_back(make_pair(10,1.));
00121   natEl.push_back(a9);
00122   vector<pair<G4int,G4double> >*b0=new vector<pair<G4int,G4double> >; // 10-Ne
00123   b0->push_back(make_pair(10,.9948));
00124   b0->push_back(make_pair(11,.9975));
00125   b0->push_back(make_pair(12,1.));
00126   natEl.push_back(b0);
00127   vector<pair<G4int,G4double> >*b1=new vector<pair<G4int,G4double> >; // 11-Na
00128   b1->push_back(make_pair(12,1.));
00129   natEl.push_back(b1);
00130   vector<pair<G4int,G4double> >*b2=new vector<pair<G4int,G4double> >; // 12-Mg
00131   b2->push_back(make_pair(12,.7899));
00132   b2->push_back(make_pair(13,.8899));
00133   b2->push_back(make_pair(14,1.));
00134   natEl.push_back(b2);
00135   vector<pair<G4int,G4double> >*b3=new vector<pair<G4int,G4double> >; // 13-Al
00136   b3->push_back(make_pair(14,1.));
00137   natEl.push_back(b3);
00138   vector<pair<G4int,G4double> >*b4=new vector<pair<G4int,G4double> >; // 14-Si
00139   b4->push_back(make_pair(14,.9223));
00140   b4->push_back(make_pair(15,.969));
00141   b4->push_back(make_pair(16,1.));
00142   natEl.push_back(b4);
00143   vector<pair<G4int,G4double> >*b5=new vector<pair<G4int,G4double> >; // 15-P
00144   b5->push_back(make_pair(16,1.));
00145   natEl.push_back(b5);
00146   vector<pair<G4int,G4double> >*b6=new vector<pair<G4int,G4double> >; // 16-S
00147   b6->push_back(make_pair(16,.9502));
00148   b6->push_back(make_pair(18,.9923));
00149   b6->push_back(make_pair(17,.9998));
00150   b6->push_back(make_pair(20,1.));
00151   natEl.push_back(b6);
00152   vector<pair<G4int,G4double> >*b7=new vector<pair<G4int,G4double> >; // 17-Cl
00153   b7->push_back(make_pair(18,.7577));
00154   b7->push_back(make_pair(20,1.));
00155   natEl.push_back(b7);
00156   vector<pair<G4int,G4double> >*b8=new vector<pair<G4int,G4double> >; // 18-Ar
00157   b8->push_back(make_pair(22,.996));
00158   b8->push_back(make_pair(18,.99937));
00159   b8->push_back(make_pair(20,1.));
00160   natEl.push_back(b8);
00161   vector<pair<G4int,G4double> >*b9=new vector<pair<G4int,G4double> >; // 19-K
00162   b9->push_back(make_pair(20,.932581));
00163   b9->push_back(make_pair(22,.999883));
00164   b9->push_back(make_pair(21,1.));
00165   natEl.push_back(b9);
00166   vector<pair<G4int,G4double> >*c0=new vector<pair<G4int,G4double> >; // 20-Ca
00167   c0->push_back(make_pair(20,.96941));
00168   c0->push_back(make_pair(24,.99027));
00169   c0->push_back(make_pair(22,.99674));
00170   c0->push_back(make_pair(28,.99861));
00171   c0->push_back(make_pair(23,.99996));
00172   c0->push_back(make_pair(26,1.));
00173   natEl.push_back(c0);
00174   vector<pair<G4int,G4double> >*c1=new vector<pair<G4int,G4double> >; // 21-Sc
00175   c1->push_back(make_pair(24,1.));
00176   natEl.push_back(c1);
00177   vector<pair<G4int,G4double> >*c2=new vector<pair<G4int,G4double> >; // 22-Ti
00178   c2->push_back(make_pair(26,.738));
00179   c2->push_back(make_pair(24,.818));
00180   c2->push_back(make_pair(25,.891));
00181   c2->push_back(make_pair(27,.946));
00182   c2->push_back(make_pair(28,1.));
00183   natEl.push_back(c2);
00184   vector<pair<G4int,G4double> >*c3=new vector<pair<G4int,G4double> >; // 23-V
00185   c3->push_back(make_pair(28,.9975));
00186   c3->push_back(make_pair(27,1.));
00187   natEl.push_back(c3);
00188   vector<pair<G4int,G4double> >*c4=new vector<pair<G4int,G4double> >; // 24-Cr
00189   c4->push_back(make_pair(28,.8379));
00190   c4->push_back(make_pair(29,.9329));
00191   c4->push_back(make_pair(26,.97635));
00192   c4->push_back(make_pair(30,1.));
00193   natEl.push_back(c4);
00194   vector<pair<G4int,G4double> >*c5=new vector<pair<G4int,G4double> >; // 25-Mn
00195   c5->push_back(make_pair(30,1.));
00196   natEl.push_back(c5);
00197   vector<pair<G4int,G4double> >*c6=new vector<pair<G4int,G4double> >; // 26-Fe
00198   c6->push_back(make_pair(30,.9172));
00199   c6->push_back(make_pair(28,.9762));
00200   c6->push_back(make_pair(31,.9972));
00201   c6->push_back(make_pair(32,1.));
00202   natEl.push_back(c6);
00203   vector<pair<G4int,G4double> >*c7=new vector<pair<G4int,G4double> >; // 27-Co
00204   c7->push_back(make_pair(32,1.));
00205   natEl.push_back(c7);
00206   vector<pair<G4int,G4double> >*c8=new vector<pair<G4int,G4double> >; // 28-Ni
00207   c8->push_back(make_pair(30,.68077));
00208   c8->push_back(make_pair(32,.943));
00209   c8->push_back(make_pair(34,.97934));
00210   c8->push_back(make_pair(33,.99074));
00211   c8->push_back(make_pair(36,1.));
00212   natEl.push_back(c8);
00213   vector<pair<G4int,G4double> >*c9=new vector<pair<G4int,G4double> >; // 29-Cu
00214   c9->push_back(make_pair(34,.6917));
00215   c9->push_back(make_pair(36,1.));
00216   natEl.push_back(c9);
00217   vector<pair<G4int,G4double> >*d0=new vector<pair<G4int,G4double> >; // 30-Zn
00218   d0->push_back(make_pair(34,.486));
00219   d0->push_back(make_pair(36,.765));
00220   d0->push_back(make_pair(38,.953));
00221   d0->push_back(make_pair(37,.994));
00222   d0->push_back(make_pair(40,1.));
00223   natEl.push_back(d0);
00224   vector<pair<G4int,G4double> >*d1=new vector<pair<G4int,G4double> >; // 31-Ga
00225   d1->push_back(make_pair(38,.60108));
00226   d1->push_back(make_pair(40,1.));
00227   natEl.push_back(d1);
00228   vector<pair<G4int,G4double> >*d2=new vector<pair<G4int,G4double> >; // 32-Ge
00229   d2->push_back(make_pair(42,.3594));
00230   d2->push_back(make_pair(40,.6360));
00231   d2->push_back(make_pair(38,.8484));
00232   d2->push_back(make_pair(41,.9256));
00233   d2->push_back(make_pair(44,1.));
00234   natEl.push_back(d2);
00235   vector<pair<G4int,G4double> >*d3=new vector<pair<G4int,G4double> >; // 33-As
00236   d3->push_back(make_pair(42,1.));
00237   natEl.push_back(d3);
00238   vector<pair<G4int,G4double> >*d4=new vector<pair<G4int,G4double> >; // 34-Se
00239   d4->push_back(make_pair(46,.4961));
00240   d4->push_back(make_pair(44,.7378));
00241   d4->push_back(make_pair(42,.8274));
00242   d4->push_back(make_pair(48,.9148));
00243   d4->push_back(make_pair(43,.9911));
00244   d4->push_back(make_pair(40,1.));
00245   natEl.push_back(d4);
00246   vector<pair<G4int,G4double> >*d5=new vector<pair<G4int,G4double> >; // 35-Br
00247   d5->push_back(make_pair(44,.5069));
00248   d5->push_back(make_pair(46,1.));
00249   natEl.push_back(d5);
00250   vector<pair<G4int,G4double> >*d6=new vector<pair<G4int,G4double> >; // 36-Kr
00251   d6->push_back(make_pair(48,.57));
00252   d6->push_back(make_pair(50,.743));
00253   d6->push_back(make_pair(46,.859));
00254   d6->push_back(make_pair(47,.974));
00255   d6->push_back(make_pair(44,.9965));
00256   d6->push_back(make_pair(42,1.));
00257   natEl.push_back(d6);
00258   vector<pair<G4int,G4double> >*d7=new vector<pair<G4int,G4double> >; // 37-Rb
00259   d7->push_back(make_pair(48,.7217));
00260   d7->push_back(make_pair(50,1.));
00261   natEl.push_back(d7);
00262   vector<pair<G4int,G4double> >*d8=new vector<pair<G4int,G4double> >; // 38-sr
00263   d8->push_back(make_pair(50,.8258));
00264   d8->push_back(make_pair(48,.9244));
00265   d8->push_back(make_pair(49,.9944));
00266   d8->push_back(make_pair(46,1.));
00267   natEl.push_back(d8);
00268   vector<pair<G4int,G4double> >*d9=new vector<pair<G4int,G4double> >; // 39-Y
00269   d9->push_back(make_pair(50,1.));
00270   natEl.push_back(d9);
00271   vector<pair<G4int,G4double> >*e0=new vector<pair<G4int,G4double> >; // 40-Zr
00272   e0->push_back(make_pair(50,.5145));
00273   e0->push_back(make_pair(54,.6883));
00274   e0->push_back(make_pair(52,.8598));
00275   e0->push_back(make_pair(51,.972));
00276   e0->push_back(make_pair(56,1.));
00277   natEl.push_back(e0);
00278   vector<pair<G4int,G4double> >*e1=new vector<pair<G4int,G4double> >; // 41-Nb
00279   e1->push_back(make_pair(52,1.));
00280   natEl.push_back(e1);
00281   vector<pair<G4int,G4double> >*e2=new vector<pair<G4int,G4double> >; // 42-Mo
00282   e2->push_back(make_pair(56,.2413));
00283   e2->push_back(make_pair(54,.4081));
00284   e2->push_back(make_pair(53,.5673));
00285   e2->push_back(make_pair(50,.7157));
00286   e2->push_back(make_pair(58,.8120));
00287   e2->push_back(make_pair(55,.9075));
00288   e2->push_back(make_pair(52,1.));
00289   natEl.push_back(e2);
00290   vector<pair<G4int,G4double> >*e3=new vector<pair<G4int,G4double> >; // 43-Tc
00291   e3->push_back(make_pair(55,1.));
00292   natEl.push_back(e3);
00293   vector<pair<G4int,G4double> >*e4=new vector<pair<G4int,G4double> >; // 44-Ru
00294   e4->push_back(make_pair(58,.316));
00295   e4->push_back(make_pair(60,.502));
00296   e4->push_back(make_pair(57,.673));
00297   e4->push_back(make_pair(55,.8));
00298   e4->push_back(make_pair(56,.926));
00299   e4->push_back(make_pair(52,.9814));
00300   e4->push_back(make_pair(54,1.));
00301   natEl.push_back(e4);
00302   vector<pair<G4int,G4double> >*e5=new vector<pair<G4int,G4double> >; // 45-Rh
00303   e5->push_back(make_pair(58,1.));
00304   natEl.push_back(e5);
00305   vector<pair<G4int,G4double> >*e6=new vector<pair<G4int,G4double> >; // 46-Pd
00306   e6->push_back(make_pair(60,.2733));
00307   e6->push_back(make_pair(62,.5379));
00308   e6->push_back(make_pair(59,.7612));
00309   e6->push_back(make_pair(55,.8784));
00310   e6->push_back(make_pair(58,.9898));
00311   e6->push_back(make_pair(56,1.));
00312   natEl.push_back(e6);
00313   vector<pair<G4int,G4double> >*e7=new vector<pair<G4int,G4double> >; // 47-Ag
00314   e7->push_back(make_pair(60,.51839));
00315   e7->push_back(make_pair(62,1.));
00316   natEl.push_back(e7);
00317   vector<pair<G4int,G4double> >*e8=new vector<pair<G4int,G4double> >; // 48-Cd
00318   e8->push_back(make_pair(66,.2873));
00319   e8->push_back(make_pair(64,.5286));
00320   e8->push_back(make_pair(59,.6566));
00321   e8->push_back(make_pair(62,.7815));
00322   e8->push_back(make_pair(65,.9037));
00323   e8->push_back(make_pair(68,.9786));
00324   e8->push_back(make_pair(58,.9911));
00325   e8->push_back(make_pair(60,1.));
00326   natEl.push_back(e8);
00327   vector<pair<G4int,G4double> >*e9=new vector<pair<G4int,G4double> >; // 49-In
00328   e9->push_back(make_pair(66,.9577));
00329   e9->push_back(make_pair(64,1.));
00330   natEl.push_back(e9);
00331   vector<pair<G4int,G4double> >*f0=new vector<pair<G4int,G4double> >; // 50-Sn
00332   f0->push_back(make_pair(70,.3259));
00333   f0->push_back(make_pair(68,.5681));
00334   f0->push_back(make_pair(66,.7134));
00335   f0->push_back(make_pair(69,.7992));
00336   f0->push_back(make_pair(67,.8760));
00337   f0->push_back(make_pair(74,.9339));
00338   f0->push_back(make_pair(72,.9802));
00339   f0->push_back(make_pair(62,.9899));
00340   f0->push_back(make_pair(64,1.));
00341   //f0->push_back(make_pair(64,.9964));
00342   //f0->push_back(make_pair(65,1.)); // Nine isotopes is the maximum, so Sn115 is out
00343   natEl.push_back(f0);
00344   vector<pair<G4int,G4double> >*f1=new vector<pair<G4int,G4double> >; // 51-Sb
00345   f1->push_back(make_pair(70,.5736));
00346   f1->push_back(make_pair(72,1.));
00347   natEl.push_back(f1);
00348   vector<pair<G4int,G4double> >*f2=new vector<pair<G4int,G4double> >; // 52-Te
00349   f2->push_back(make_pair(78,.3387));
00350   f2->push_back(make_pair(76,.6557));
00351   f2->push_back(make_pair(74,.8450));
00352   f2->push_back(make_pair(73,.9162));
00353   f2->push_back(make_pair(72,.9641));
00354   f2->push_back(make_pair(70,.9900));
00355   f2->push_back(make_pair(71,.99905));
00356   f2->push_back(make_pair(68,1.));
00357   natEl.push_back(f2);
00358   vector<pair<G4int,G4double> >*f3=new vector<pair<G4int,G4double> >; // 53-I
00359   f3->push_back(make_pair(74,1.));
00360   natEl.push_back(f3);
00361   vector<pair<G4int,G4double> >*f4=new vector<pair<G4int,G4double> >; // 54-Xe
00362   f4->push_back(make_pair(78,.269));
00363   f4->push_back(make_pair(75,.533));
00364   f4->push_back(make_pair(77,.745));
00365   f4->push_back(make_pair(80,.849));
00366   f4->push_back(make_pair(82,.938));
00367   f4->push_back(make_pair(76,.979));
00368   f4->push_back(make_pair(74,.9981));
00369   f4->push_back(make_pair(70,.9991));
00370   f4->push_back(make_pair(72,1.));
00371   natEl.push_back(f4);
00372   vector<pair<G4int,G4double> >*f5=new vector<pair<G4int,G4double> >; // 55-Cs
00373   f5->push_back(make_pair(78,1.));
00374   natEl.push_back(f5);
00375   vector<pair<G4int,G4double> >*f6=new vector<pair<G4int,G4double> >; // 56-Ba
00376   f6->push_back(make_pair(82,.717));
00377   f6->push_back(make_pair(81,.8293));
00378   f6->push_back(make_pair(80,.9078));
00379   f6->push_back(make_pair(79,.97373));
00380   f6->push_back(make_pair(78,.99793));
00381   f6->push_back(make_pair(74,.99899));
00382   f6->push_back(make_pair(76,1.));
00383   natEl.push_back(f6);
00384   vector<pair<G4int,G4double> >*f7=new vector<pair<G4int,G4double> >; // 57-La
00385   f7->push_back(make_pair(82,.999098));
00386   f7->push_back(make_pair(81,1.));
00387   natEl.push_back(f7);
00388   vector<pair<G4int,G4double> >*f8=new vector<pair<G4int,G4double> >; // 58-Ce
00389   f8->push_back(make_pair(82,.8843));
00390   f8->push_back(make_pair(84,.9956));
00391   f8->push_back(make_pair(80,.9981));
00392   f8->push_back(make_pair(78,1.));
00393   natEl.push_back(f8);
00394   vector<pair<G4int,G4double> >*f9=new vector<pair<G4int,G4double> >; // 59-Pr
00395   f9->push_back(make_pair(82,1.));
00396   natEl.push_back(f9);
00397   vector<pair<G4int,G4double> >*g0=new vector<pair<G4int,G4double> >; // 60-Nd
00398   g0->push_back(make_pair(82,.2713));
00399   g0->push_back(make_pair(84,.5093));
00400   g0->push_back(make_pair(86,.6812));
00401   g0->push_back(make_pair(83,.8030));
00402   g0->push_back(make_pair(85,.8860));
00403   g0->push_back(make_pair(88,.9436));
00404   g0->push_back(make_pair(90,1.));
00405   natEl.push_back(g0);
00406   vector<pair<G4int,G4double> >*g1=new vector<pair<G4int,G4double> >; // 61-Pm
00407   g1->push_back(make_pair(85,1.));
00408   natEl.push_back(g1);
00409   vector<pair<G4int,G4double> >*g2=new vector<pair<G4int,G4double> >; // 62-Sm
00410   g2->push_back(make_pair(90,.267));
00411   g2->push_back(make_pair(92,.494));
00412   g2->push_back(make_pair(85,.644));
00413   g2->push_back(make_pair(87,.782));
00414   g2->push_back(make_pair(86,.895));
00415   g2->push_back(make_pair(88,.969));
00416   g2->push_back(make_pair(82,1.));
00417   natEl.push_back(g2);
00418   vector<pair<G4int,G4double> >*g3=new vector<pair<G4int,G4double> >; // 63-Eu
00419   g3->push_back(make_pair(90,.522));
00420   g3->push_back(make_pair(89,1.));
00421   natEl.push_back(g3);
00422   vector<pair<G4int,G4double> >*g4=new vector<pair<G4int,G4double> >; // 64-Gd
00423   g4->push_back(make_pair(94,.2484));
00424   g4->push_back(make_pair(96,.4670));
00425   g4->push_back(make_pair(92,.6717));
00426   g4->push_back(make_pair(93,.8282));
00427   g4->push_back(make_pair(91,.9762));
00428   g4->push_back(make_pair(90,.9980));
00429   g4->push_back(make_pair(88,1.));
00430   natEl.push_back(g4);
00431   vector<pair<G4int,G4double> >*g5=new vector<pair<G4int,G4double> >; // 65-Tb
00432   g5->push_back(make_pair(94,1.));
00433   natEl.push_back(g5);
00434   vector<pair<G4int,G4double> >*g6=new vector<pair<G4int,G4double> >; // 66-Dy
00435   g6->push_back(make_pair(98,.282));
00436   g6->push_back(make_pair(96,.537));
00437   g6->push_back(make_pair(97,.786));
00438   g6->push_back(make_pair(95,.975));
00439   g6->push_back(make_pair(94,.9984));
00440   g6->push_back(make_pair(92,.9994));
00441   g6->push_back(make_pair(90,1.));
00442   natEl.push_back(g6);
00443   vector<pair<G4int,G4double> >*g7=new vector<pair<G4int,G4double> >; // 67-Ho
00444   g7->push_back(make_pair(98,1.));
00445   natEl.push_back(g7);
00446   vector<pair<G4int,G4double> >*g8=new vector<pair<G4int,G4double> >; // 68-Er
00447   g8->push_back(make_pair( 98,.3360));
00448   g8->push_back(make_pair(100,.6040));
00449   g8->push_back(make_pair( 99,.8335));
00450   g8->push_back(make_pair(102,.9825));
00451   g8->push_back(make_pair( 96,.9986));
00452   g8->push_back(make_pair( 94,1.));
00453   natEl.push_back(g8);
00454   vector<pair<G4int,G4double> >*g9=new vector<pair<G4int,G4double> >; // 69-Tm
00455   g9->push_back(make_pair(100,1.));
00456   natEl.push_back(g9);
00457   vector<pair<G4int,G4double> >*h0=new vector<pair<G4int,G4double> >; // 70-Yb
00458   h0->push_back(make_pair(104,.3180));
00459   h0->push_back(make_pair(102,.5370));
00460   h0->push_back(make_pair(103,.6982));
00461   h0->push_back(make_pair(101,.8412));
00462   h0->push_back(make_pair(106,.9682));
00463   h0->push_back(make_pair(100,.9987));
00464   h0->push_back(make_pair( 98,1.));
00465   natEl.push_back(h0);
00466   vector<pair<G4int,G4double> >*h1=new vector<pair<G4int,G4double> >; // 71-Lu
00467   h1->push_back(make_pair(104,.9741));
00468   h1->push_back(make_pair(105,1.));
00469   natEl.push_back(h1);
00470   vector<pair<G4int,G4double> >*h2=new vector<pair<G4int,G4double> >; // 72-Hf
00471   h2->push_back(make_pair(108,.35100));
00472   h2->push_back(make_pair(106,.62397));
00473   h2->push_back(make_pair(105,.81003));
00474   h2->push_back(make_pair(107,.94632));
00475   h2->push_back(make_pair(104,.99838));
00476   h2->push_back(make_pair(102,1.));
00477   natEl.push_back(h2);
00478   vector<pair<G4int,G4double> >*h3=new vector<pair<G4int,G4double> >; // 73-Ta
00479   h3->push_back(make_pair(108,.99988));
00480   h3->push_back(make_pair(107,1.));
00481   natEl.push_back(h3);
00482   vector<pair<G4int,G4double> >*h4=new vector<pair<G4int,G4double> >; // 74-W
00483   h4->push_back(make_pair(110,.307));
00484   h4->push_back(make_pair(112,.593));
00485   h4->push_back(make_pair(108,.856));
00486   h4->push_back(make_pair(109,.9988));
00487   h4->push_back(make_pair(106,1.));
00488   natEl.push_back(h4);
00489   vector<pair<G4int,G4double> >*h5=new vector<pair<G4int,G4double> >; // 75-Re
00490   h5->push_back(make_pair(112,.626));
00491   h5->push_back(make_pair(110,1.));
00492   natEl.push_back(h5);
00493   vector<pair<G4int,G4double> >*h6=new vector<pair<G4int,G4double> >; // 78-Os
00494   h6->push_back(make_pair(116,.410));
00495   h6->push_back(make_pair(114,.674));
00496   h6->push_back(make_pair(113,.835));
00497   h6->push_back(make_pair(112,.968));
00498   h6->push_back(make_pair(111,.984));
00499   h6->push_back(make_pair(110,.9998));
00500   h6->push_back(make_pair(108,1.));
00501   natEl.push_back(h6);
00502   vector<pair<G4int,G4double> >*h7=new vector<pair<G4int,G4double> >; // 77-Ir
00503   h7->push_back(make_pair(116,.627));
00504   h7->push_back(make_pair(114,1.));
00505   natEl.push_back(h7);
00506   vector<pair<G4int,G4double> >*h8=new vector<pair<G4int,G4double> >; // 78-Pt
00507   h8->push_back(make_pair(117,.338));
00508   h8->push_back(make_pair(116,.667));
00509   h8->push_back(make_pair(118,.920));
00510   h8->push_back(make_pair(120,.992));
00511   h8->push_back(make_pair(114,.9999));
00512   h8->push_back(make_pair(112,1.));
00513   natEl.push_back(h8);
00514   vector<pair<G4int,G4double> >*h9=new vector<pair<G4int,G4double> >; // 79-Au
00515   h9->push_back(make_pair(118,1.));
00516   natEl.push_back(h9);
00517   vector<pair<G4int,G4double> >*i0=new vector<pair<G4int,G4double> >; // 80-Hg
00518   i0->push_back(make_pair(122,.2986));
00519   i0->push_back(make_pair(120,.5296));
00520   i0->push_back(make_pair(119,.6983));
00521   i0->push_back(make_pair(121,.8301));
00522   i0->push_back(make_pair(118,.9298));
00523   i0->push_back(make_pair(124,.9985));
00524   i0->push_back(make_pair(116,1.));
00525   natEl.push_back(i0);
00526   vector<pair<G4int,G4double> >*i1=new vector<pair<G4int,G4double> >; // 81-Tl
00527   i1->push_back(make_pair(124,.70476));
00528   i1->push_back(make_pair(122,1.));
00529   natEl.push_back(i1);
00530   vector<pair<G4int,G4double> >*i2=new vector<pair<G4int,G4double> >; // 82-Pb
00531   i2->push_back(make_pair(126,.524));
00532   i2->push_back(make_pair(124,.765));
00533   i2->push_back(make_pair(125,.986));
00534   i2->push_back(make_pair(122,1.));
00535   natEl.push_back(i2);
00536   vector<pair<G4int,G4double> >*i3=new vector<pair<G4int,G4double> >; // 83-Bi
00537   i3->push_back(make_pair(126,1.));
00538   natEl.push_back(i3);
00539   vector<pair<G4int,G4double> >*i4=new vector<pair<G4int,G4double> >; // 84-Po
00540   i4->push_back(make_pair(125,1.));
00541   natEl.push_back(i4);
00542   vector<pair<G4int,G4double> >*i5=new vector<pair<G4int,G4double> >; // 85-At
00543   i5->push_back(make_pair(136,1.));
00544   natEl.push_back(i5);
00545   vector<pair<G4int,G4double> >*i6=new vector<pair<G4int,G4double> >; // 86-Ru
00546   i6->push_back(make_pair(136,1.));
00547   natEl.push_back(i6);
00548   vector<pair<G4int,G4double> >*i7=new vector<pair<G4int,G4double> >; // 87-Fr
00549   i7->push_back(make_pair(138,1.));
00550   natEl.push_back(i7);
00551   vector<pair<G4int,G4double> >*i8=new vector<pair<G4int,G4double> >; // 88-Ra
00552   i8->push_back(make_pair(138,1.));
00553   natEl.push_back(i8);
00554   vector<pair<G4int,G4double> >*i9=new vector<pair<G4int,G4double> >; // 89-Ac
00555   i9->push_back(make_pair(142,1.));
00556   natEl.push_back(i9);
00557   vector<pair<G4int,G4double> >*j0=new vector<pair<G4int,G4double> >; // 90-Th
00558   j0->push_back(make_pair(142,1.));
00559   natEl.push_back(j0);
00560   vector<pair<G4int,G4double> >*j1=new vector<pair<G4int,G4double> >; // 91-Pa
00561   j1->push_back(make_pair(140,1.));
00562   natEl.push_back(j1);
00563   vector<pair<G4int,G4double> >*j2=new vector<pair<G4int,G4double> >; // 92-U
00564   j2->push_back(make_pair(146,.992745));
00565   j2->push_back(make_pair(143,.999945));
00566   j2->push_back(make_pair(142,1.));
00567   natEl.push_back(j2);
00568   vector<pair<G4int,G4double> >*j3=new vector<pair<G4int,G4double> >; // 93-Np
00569   j3->push_back(make_pair(144,1.));
00570   natEl.push_back(j3);
00571   vector<pair<G4int,G4double> >*j4=new vector<pair<G4int,G4double> >; // 94-Pu
00572   j4->push_back(make_pair(150,1.));
00573   natEl.push_back(j4);
00574   vector<pair<G4int,G4double> >*j5=new vector<pair<G4int,G4double> >; // 95-Am
00575   j5->push_back(make_pair(148,1.));
00576   natEl.push_back(j5);
00577   vector<pair<G4int,G4double> >*j6=new vector<pair<G4int,G4double> >; // 96-Cm
00578   j6->push_back(make_pair(151,1.));
00579   natEl.push_back(j6);
00580   vector<pair<G4int,G4double> >*j7=new vector<pair<G4int,G4double> >; // 97-Bk
00581   j7->push_back(make_pair(150,1.));
00582   natEl.push_back(j7);
00583   vector<pair<G4int,G4double> >*j8=new vector<pair<G4int,G4double> >; // 98-Cf
00584   j8->push_back(make_pair(153,1.));
00585   natEl.push_back(j8);
00586   vector<pair<G4int,G4double> >*j9=new vector<pair<G4int,G4double> >; // 99-Es
00587   j9->push_back(make_pair(157,1.));
00588   natEl.push_back(j9);
00589   vector<pair<G4int,G4double> >*k0=new vector<pair<G4int,G4double> >; // 100-Fm
00590   k0->push_back(make_pair(157,1.));
00591   natEl.push_back(k0);
00592   vector<pair<G4int,G4double> >*k1=new vector<pair<G4int,G4double> >; // 101-Md
00593   k1->push_back(make_pair(157,1.));
00594   natEl.push_back(k1);
00595   vector<pair<G4int,G4double> >*k2=new vector<pair<G4int,G4double> >; // 102-No
00596   k2->push_back(make_pair(157,1.));
00597   natEl.push_back(k2);
00598   vector<pair<G4int,G4double> >*k3=new vector<pair<G4int,G4double> >; // 103-Lr
00599   k3->push_back(make_pair(157,1.));
00600   natEl.push_back(k3);
00601   vector<pair<G4int,G4double> >*k4=new vector<pair<G4int,G4double> >; // 104-Rf
00602   k4->push_back(make_pair(157,1.));
00603   natEl.push_back(k4);
00604   vector<pair<G4int,G4double> >*k5=new vector<pair<G4int,G4double> >; // 105-Db
00605   k5->push_back(make_pair(157,1.));
00606   natEl.push_back(k5);
00607   vector<pair<G4int,G4double> >*k6=new vector<pair<G4int,G4double> >; // 106-Sg
00608   k6->push_back(make_pair(157,1.));
00609   natEl.push_back(k6);
00610   vector<pair<G4int,G4double> >*k7=new vector<pair<G4int,G4double> >; // 107-Bh
00611   k7->push_back(make_pair(155,1.));
00612   natEl.push_back(k7);
00613   vector<pair<G4int,G4double> >*k8=new vector<pair<G4int,G4double> >; // 108-Hs
00614   k8->push_back(make_pair(157,1.));
00615   natEl.push_back(k8);
00616   vector<pair<G4int,G4double> >*k9=new vector<pair<G4int,G4double> >; // 109-Mt
00617   k9->push_back(make_pair(157,1.));
00618   natEl.push_back(k9);
00619   // Now fill natElements and natIsoCrossS
00620   G4int nona=natEl.size();
00621 #ifdef cdebug
00622   G4cout<<"G4QIsotope::Constructor natEl filling is finished nE="<<nona<<G4endl;
00623 #endif
00624   for(G4int i=0; i<nona; i++)
00625   {
00626     vector<pair<G4int,G4double> >* is=natEl[i]; // Pointer to theElement
00627     G4int n=is->size();
00628 #ifdef cdebug
00629     G4cout<<"G4QIsotope::Constructor: Element # "<<i<<", nOfIsotopes="<<n<<G4endl;
00630 #endif
00631     vector<pair<G4int,G4double>*>*a=new vector<pair<G4int,G4double>*>;
00632     vector<pair<G4int,G4double>*>*s_vec=new vector<pair<G4int,G4double>*>;
00633     G4double last=0.;
00634     if(n) for(G4int j=0; j<n; j++)
00635     {
00636       G4int    nn =(*is)[j].first; // #ofNeutrons in the isotope
00637       G4double cur=(*is)[j].second;// value of the summed abundancy
00638       pair<G4int,G4double>* aP = new pair<G4int,G4double>(nn,cur-last);
00639       last=cur;                     // Update the summed value
00640       pair<G4int,G4double>* sP = new pair<G4int,G4double>((*is)[j]);
00641       a->push_back(aP);
00642       s_vec->push_back(sP);
00643 #ifdef cdebug
00644       G4cout<<"G4QIsotope::Constructor:Element# "<<i<<", Pair # "<<j<<" is filled"<<G4endl;
00645 #endif
00646     }
00647     natElements.push_back(a);       // Fill abundancies for the particular isotope
00648     natSumAbund.push_back(s_vec);   // Fill summes abundancies up to this isotope
00649 #ifdef cdebug
00650     G4cout<<"G4QIsotope::Constructor: natElements is filled"<<G4endl;
00651 #endif
00652     vector<pair<G4int,G4double>*>*c=new vector<pair<G4int,G4double>*>;
00653     if(n) for(G4int j=0; j<n; j++)  // Cross sections are 0. by default
00654     {
00655       pair<G4int,G4double>* cP = new pair<G4int,G4double>((*is)[j].first,0.);
00656       c->push_back(cP);
00657 #ifdef cdebug
00658       G4cout<<"G4QIsotope::Constructor:CrosSecPair i="<<i<<", j="<<j<<" is filled"<<G4endl;
00659 #endif
00660     }
00661     natIsoCrosS.push_back(c);       // FillPrototypeCrossSec's (0) for theParticularIsotope
00662 #ifdef cdebug
00663     G4cout<<"G4QIsotope::Constructor: natIsoCrosS is filled"<<G4endl;
00664 #endif
00665     delete is;
00666   }
00667 #ifdef cdebug
00668   G4cout<<"G4QIsotope::Constructor: is finished"<<G4endl;
00669 #endif
00670 }
00671 
00672 G4QIsotope::~G4QIsotope()          // The QIsotopes are destructed only in theEnd of theJob
00673 {
00674 #ifdef debug
00675   G4cout<<"G4QIsotope::Destructor is called"<<G4endl;
00676 #endif
00677   G4int uP=natElements.size();
00678   if(uP) for(G4int i=0; i<uP; i++)
00679   {
00680     vector<pair<G4int,G4double>*>* curA=natElements[i];
00681     G4int nn=curA->size();         // Can not be 0 by definition
00682     if(nn) for(G4int n=0; n<nn; n++) delete (*curA)[n]; // Delete pair(N,Ab)
00683     delete curA;                   // Delet abundancy vector
00684     vector<pair<G4int,G4double>*>* curS=natSumAbund[i];
00685     G4int ns_value=curS->size();   // Can not be 0 by definition
00686     if(ns_value) for(G4int n=0; n<ns_value; n++) delete (*curS)[n]; // Delete pair(N,Ab)
00687     delete curS;                   // Delet abundancy vector
00688     vector<pair<G4int,G4double>*>* curC=natIsoCrosS[i];
00689     G4int nc=curC->size();         // Can not be 0 by definition
00690     if(nc) for(G4int k=0; k<nc; k++) delete (*curC)[k]; // Delete pair(N,CS)
00691     delete curC;                   // Delete cross section vector
00692   }
00693   G4int nP=newElems.size();
00694   if(nP) for(G4int j=0; j<nP; j++) // LOOP over new UserDefinedElements
00695   {
00696     pair<G4int, vector<pair<G4int,G4double>*>* >* nEl= newElems[j];
00697     G4int nEn=nEl->second->size();
00698     if(nEn) for(G4int k=0; k<nEn; k++) delete (*(nEl->second))[k]; // Del vect<pair(N,A)*>
00699     delete nEl->second;            // Delete the vector
00700     delete nEl;                    // Delete vect<IndZ,vect<pair(N,Ab)*>*> newElementVector
00701     //
00702     pair<G4int, vector<pair<G4int,G4double>*>* >* nSA= newSumAb[j];
00703     G4int nSn=nSA->second->size();
00704     if(nSn) for(G4int n=0; n<nSn; n++) delete (*(nSA->second))[n]; // Del vect<pair(N,S)*>
00705     delete nSA->second;            // Delete the vector
00706     delete nSA;                    // Delete vect<IndZ,vect<pair(N,SA)*>*> newSumAbunVector
00707     //
00708     pair<G4int, vector<pair<G4int,G4double>*>* >* nCS= newIsoCS[j];
00709     G4int nCn=nCS->second->size();
00710     if(nCn) for(G4int n=0; n<nCn; n++) delete (*(nCS->second))[n]; // Del vect<pair(N,C)*>
00711     delete nCS->second;            // Delete the vector
00712     delete nCS;                    // Delete vect<IndZ,vect<pair(N,CS)*>*> newIsoCroSVector
00713     //
00714     if(nEn!=nCn) G4cerr<<"*G4QIsotope-WORNING-:#El="<<j<<":nE="<<nEn<<"!=nC="<<nCn<<G4endl;
00715     if(nEn!=nSn) G4cerr<<"*G4QIsotope-WORNING-:#El="<<j<<":nE="<<nEn<<"!=nS="<<nSn<<G4endl;
00716   }
00717 }
00718 
00719 // Returns Pointer to the G4QIsotope singletone
00720 G4QIsotope* G4QIsotope::Get()
00721 {
00722 #ifdef pdebug
00723   G4cout<<"G4QIsotope::Get is called"<<G4endl;
00724 #endif
00725   static G4QIsotope theIsotopes;             // *** Static body of the G4QIsotope class ***
00726   return &theIsotopes;
00727 }
00728 
00729 // #ofProtons in stable isotopes with fixed A=Z+N. Returns length and fils VectOfIsotopes
00730 G4int G4QIsotope::GetProtons(G4int A, vector<G4int>& isoV) 
00731 {
00732   const G4int nAZ=270;  // Dimension of the table
00733   // Best Z for the given A
00734   const G4int bestZ[nAZ] = {
00735      0,  1,  1,  2,  2,  0,  3,  3,  4,  4,   //0
00736      5,  5,  6,  6,  7,  7,  8,  8,  8,  9,   //10
00737     10, 10, 10, 11, 12, 12, 12, 13, 14, 14,   //20
00738     14, 15, 16, 16, 16, 17, 18, 17, 18, 19,   //30
00739     18, 19, 20, 20, 20, 21, 22, 22, 23, 23,   //40
00740     22, 23, 24, 24, 26, 25, 26, 26, 28, 27,   //50
00741     28, 28, 28, 29, 30, 29, 30, 30, 30, 31,   //60
00742     32, 31, 32, 32, 32, 33, 34, 34, 34, 35,   //70
00743     34, 35, 36, 36, 36, 37, 39, 36, 38, 39,   //80
00744     40, 40, 41, 40, 40, 42, 42, 42, 42, 44,   //90
00745     44, 44, 44, 45, 44, 46, 46, 47, 46, 47,   //100
00746     48, 48, 48, 48, 48, 49, 50, 50, 50, 50,   //110
00747     50, 51, 50, 51, 50, 52, 52, 53, 52, 54,   //120
00748     52, 54, 54, 55, 54, 56, 54, 56, 56, 57,   //130
00749     58, 59, 60, 60, 60, 60, 60, 62, 62, 62,   //140
00750     62, 63, 62, 63, 62, 64, 64, 64, 64, 65,   //150
00751     64, 66, 66, 66, 66, 67, 68, 68, 68, 69,   //160
00752     68, 70, 70, 70, 70, 71, 70, 72, 72, 72,   //170
00753     72, 73, 74, 74, 74, 75, 74, 75, 76, 76,   //180
00754     76, 77, 76, 77, 78, 78, 78, 79, 80, 80,   //190
00755     80, 80, 80, 81, 80, 81, 82, 82, 82, 83,   //200
00756     82,  0, 82,  0, 82,  0, 84,  0,  0,  0,   //210
00757     86,  0, 86, 87, 88,  0, 88, 89, 88, 89,   //220
00758     89, 91, 90,  0, 92, 92,  0, 93, 92, 94,   //230
00759      0,  0,  0, 95, 94,  0,  0, 96,  0,  0,   //240
00760      0, 98, 99,  0,  0,  0,  0,100,101,102,   //250
00761    103,104,105,106,  0,108,109,  0,  0,  0};  //260
00762   // 0   1   2   3   4   5   6   7   8   9
00763   // Second candidate
00764   const G4int secoZ[nAZ] = {
00765      0,  0,  0,  1,  0,  0,  0,  4,  0,  0,   //0
00766      4,  6,  5,  7,  0,  8,  0,  0,  0,  8,   //10
00767      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,   //20
00768      0,  0, 15,  0,  0, 16,  0,  0,  0,  0,   //30
00769     20, 20,  0,  0,  0,  0, 20,  0,  0,  0,   //40
00770     24,  0,  0,  0, 24,  0,  0,  0, 26,  0,   //50
00771     27,  0,  0,  0, 28,  0,  0,  0,  0,  0,   //60
00772     30,  0,  0,  0, 34,  0, 32,  0,  0,  0,   //70
00773     36,  0, 34,  0, 38,  0, 38, 38,  0,  0,   //80
00774      0,  0, 42,  0, 42,  0, 44,  0, 44,  0,   //90
00775     42,  0, 46,  0, 46,  0, 48,  0, 48,  0,   //100
00776     46,  0, 50, 49, 50, 50, 48,  0,  0,  0,   //110
00777     52,  0, 52,  0, 52,  0, 54,  0, 54,  0,   //120
00778     54,  0, 56,  0, 56,  0, 56,  0, 58,  0,   //130
00779     54,  0, 58,  0, 62, 61,  0,  0, 60,  0,   //140
00780     60,  0, 64,  0, 64,  0, 66,  0, 66,  0,   //150
00781     66,  0, 68,  0, 68,  0,  0,  0, 70,  0,   //160
00782     70,  0,  0,  0, 72,  0, 72,  0,  0,  0,   //170
00783     74,  0,  0,  0, 76,  0, 76, 76,  0,  0,   //180
00784     78,  0, 78,  0,  0,  0, 80,  0, 78,  0,   //190
00785      0,  0,  0,  0, 82,  0,  0,  0,  0, 84,   //200
00786     84,  0, 83,  0, 83,  0,  0,  0,  0,  0,   //210
00787      0,  0,  0,  0,  0,  0,  0,  0, 89,  0,   //220
00788      0,  0,  0,  0, 93,  0,  0,  0, 93,  0,   //230
00789      0,  0,  0,  0,  0,  0,  0, 97,  0,  0,   //240
00790      0,  0,  0,  0,  0,  0,  0,  0,  0,  0,   //250
00791      0,  0,107,  0,  0,  0,  0,  0,  0,  0};  //260
00792   // 0   1   2   3   4   5   6   7   8   9
00793   // Third candidate
00794   const G4int thrdZ[nAZ] = {
00795     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //0
00796     0, 0, 7, 0, 0, 0, 0, 0, 0, 0,   //10
00797     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //20
00798     0, 0, 0, 0, 0,20, 0, 0, 0, 0,   //30
00799    19, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //40
00800    23, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //50
00801     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //60
00802     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //70
00803     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //80
00804     0, 0,36, 0,38, 0,40, 0, 0, 0,   //90
00805     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //100
00806     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //110
00807     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //120
00808    56, 0,50, 0, 0, 0,58, 0,57, 0,   //130
00809     0, 0, 0, 0, 0, 0, 0, 0,65, 0,   //140
00810     0, 0,66, 0, 0, 0, 0, 0, 0, 0,   //150
00811     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //160
00812     0, 0, 0, 0, 0, 0,71, 0, 0, 0,   //170
00813    73, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //180
00814     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //190
00815     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //200
00816    83, 0,84, 0,84, 0, 0, 0, 0, 0,   //210
00817     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //220
00818     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //230
00819     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //240
00820     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //250
00821     0, 0, 0, 0, 0, 0, 0, 0, 0, 0};  //260
00822   //0  1  2  3  4  5  6  7  8  9
00823   // Fourth candidate (only two isotopes)
00824   const G4int quadZ[nAZ] = {
00825     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //0
00826     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //10
00827     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //20
00828     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //30
00829     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //40
00830     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //50
00831     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //60
00832     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //70
00833     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //80
00834     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //90
00835     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //100
00836     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //110
00837     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //120
00838     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //130
00839     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //140
00840     0, 0,67, 0, 0, 0, 0, 0, 0, 0,   //150
00841     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //160
00842     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //170
00843     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //180
00844     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //190
00845     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //200
00846    85, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //210
00847     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //220
00848     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //230
00849     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //240
00850     0, 0, 0, 0, 0, 0, 0, 0, 0, 0,   //250
00851     0, 0, 0, 0, 0, 0, 0, 0, 0, 0};  //260
00852   //0  1  2  3  4  5  6  7  8  9
00853   isoV.clear();                     // Always cleans up before filling
00854   if(A>=nAZ) return 0;
00855   G4int fZ=bestZ[A];
00856   if(fZ)
00857   {
00858     isoV.push_back(fZ);
00859     G4int sZ=secoZ[A];
00860     if(sZ)
00861     {
00862       isoV.push_back(sZ);
00863       G4int tZ=thrdZ[A];
00864       if(tZ)
00865       {
00866         isoV.push_back(tZ);
00867         G4int qZ=quadZ[A];
00868         if(qZ)
00869         {
00870           isoV.push_back(qZ);
00871           return 4;
00872         }
00873         else return 3;
00874       }
00875       else return 2;
00876     }
00877     else return 1;
00878   }
00879   else return 0;
00880 }
00881 
00882 // Randomize a#ofNeutrons in the Element (independent function @@ can be used for initial)
00883 G4int G4QIsotope::RandomizeNeutrons(G4int i) 
00884 {
00885   static const G4int nElements = 110; // Max=Meitnerium(Mt)Z=99(starts with Z=0 - neuteron)
00886   // If an error is found in this simple tree, please, correct the initialization above
00887   static G4int dN[nElements]=
00888   {
00889  //   n              Be                   F              Al       P
00890       1, -1, -1, -1,  5, -1, -1, -1, -1, 10, -1, 12, -1, 14, -1, 16, -1, -1, -1, -1,
00891  //      Sc              Mn      Co                      As                       Y
00892      -1, 24, -1, -1, -1, 30, -1, 32, -1, -1, -1, -1, -1, 42, -1, -1, -1, -1, -1, 50,
00893  //      Nb      Tc      Rh                               I      Cs              Pr
00894      -1, 52, -1, 55, -1, 58, -1, -1, -1, -1, -1, -1, -1, 74, -1, 78, -1, -1, -1, 82,
00895  //      Pm              Tb      Ho      Tm                                      Au
00896      -1,-85, -1, -1, -1, 94, -1, 98, -1,100, -1, -1, -1, -1, -1, -1, -1, -1, -1,118,
00897  //                      Po   At   Rn   Fr   Ra   Ac
00898      -1,  -1,  -1, 126,-125,-136,-136,-138,-138,-142,
00899  //  Th   Pa    U   Np   Pu   Am   Cm   Bk   Cf   Es
00900     142,-140,  -1,-144,-150,-148,-151,-150,-153,-153
00901  //  Fm   Md   No   Lr   Rf   Db   Sg   Bh   Hs   Mt
00902    -157,-157,-157,-157,-157,-157,-157,-155,-157,-157
00903   };
00904 #ifdef debug
00905   G4cout<<"G4QIsotope::RandomizeNeutrons is called Z="<<i<<G4endl;
00906 #endif
00907   G4int N=dN[i];
00908   if(N==-1)
00909   {
00910     G4double rnd=G4UniformRand();
00911     if          (i<44)        // =----= H - Mo
00912     {
00913       if        (i<23)        // ------ H - Ti
00914       {
00915         if      (i<12)        // ______ H - Ne
00916         {
00917           if     (i<6)        // ...... H - B
00918           {
00919             if   (i<3)
00920             {
00921               if(i==1)        // H
00922               {
00923                 if(rnd>.00015)       N=0;
00924                 else                 N=1;
00925               }
00926               else            // He (2)
00927               {
00928                 if(rnd>1.37e-6)      N=2;
00929                 else                 N=1;
00930               }
00931             }
00932             else
00933             {
00934               if(i==3)        // Li
00935               {
00936                 if(rnd>.075)         N=4;
00937                 else                 N=3;
00938               }
00939               else            // B (5)
00940               {
00941                 if(rnd>.199)         N=6;
00942                 else                 N=5;
00943               }
00944             }
00945           }
00946           else                // ...... C - Ne
00947           {
00948             if   (i<8)
00949             {
00950               if(i==6)        // C
00951               {
00952                 if(rnd>.011)         N=6;
00953                 else                 N=7;
00954               }
00955               else            // N (7)
00956               {
00957                 if(rnd>.0037)        N=7;
00958                 else                 N=8;
00959               }
00960             }
00961             else
00962             {
00963               if(i==8)        // O
00964               {
00965                 if     (rnd<.9976)   N=8;
00966                 else if(rnd<.9996)   N=10;
00967                 else                 N=9;
00968               }
00969               else            // Ne (10)
00970               {
00971                 if     (rnd<.9948)   N=10;
00972                 else if(rnd<.9975)   N=11;
00973                 else                 N=12;
00974               }
00975             }
00976           }
00977         }
00978         else                  // ______ Mg - Ti
00979         {
00980           if     (i<18)       // ...... Mg - Cl
00981           {
00982             if   (i<16)
00983             {
00984               if(i==12)       // Mg
00985               {
00986                 if     (rnd<.7899)   N=12;
00987                 else if(rnd<.8899)   N=13;
00988                 else                 N=14;
00989               }
00990               else            // Si (14)
00991               {
00992                 if    (rnd<.9223)    N=14;
00993                 else if(rnd<.969)    N=15;
00994                 else                 N=16;
00995               }
00996             }
00997             else
00998             {
00999               if(i==16)       // S
01000               {
01001                 if     (rnd<.9502)   N=16;
01002                 else if(rnd<.9923)   N=18;
01003                 else if(rnd<.9998)   N=17;
01004                 else                 N=20;
01005               }
01006               else            // Cl (17)
01007               {
01008                 if     (rnd>.7577)   N=18;
01009                 else                 N=20;
01010               }
01011             }
01012           }
01013           else                // ...... Ar - Ti
01014           {
01015             if   (i<20)
01016             {
01017               if(i==18)       // Ar
01018               {
01019                 if     (rnd<.996)    N=22;
01020                 else if(rnd<.99937)  N=18;
01021                 else                 N=20;
01022               }
01023               else            // K (19)
01024               {
01025                 if     (rnd<.932581) N=20;
01026                 else if(rnd<.999883) N=22;
01027                 else                 N=21;
01028               }
01029             }
01030             else
01031             {
01032               if(i==20)       // Ca
01033               {
01034                 if     (rnd<.96941)  N=20;
01035                 else if(rnd<.99027)  N=24;
01036                 else if(rnd<.99674)  N=22;
01037                 else if(rnd<.99861)  N=28;
01038                 else if(rnd<.99996)  N=23;
01039                 else                 N=26;
01040               }
01041               else            // Ti (22)
01042               {
01043                 if     (rnd<.738)    N=26;
01044                 else if(rnd<.818)    N=24;
01045                 else if(rnd<.891)    N=25;
01046                 else if(rnd<.946)    N=27;
01047                 else                 N=28;
01048               }
01049             }
01050           }
01051         }
01052       }
01053       else                    // ------ V - Mo
01054       {
01055         if     (i<32)         // ______ V - Ga
01056         {
01057           if     (i<28)       // ...... H - Fe
01058           {
01059             if   (i<26)
01060             {
01061               if(i==23)       // V
01062               {
01063                 if     (rnd<.9975)   N=28;
01064                 else                 N=27;
01065               }
01066               else            // Cr (24)
01067               {
01068                 if     (rnd<.8379)   N=28;
01069                 else if(rnd<.9329)   N=29;
01070                 else if(rnd<.97635)  N=26;
01071                 else                 N=30;
01072               }
01073             }
01074             else              // Fe (26)
01075             {
01076                 if     (rnd<.9172)   N=30;
01077                 else if(rnd<.9762)   N=28;
01078                 else if(rnd<.9972)   N=31;
01079                 else                 N=32;
01080             }
01081           }
01082           else                // ...... Ni - Ga
01083           {
01084             if   (i<30)
01085             {
01086               if(i==28)       // Ni
01087               {
01088                 if     (rnd<.68077)  N=30;
01089                 else if(rnd<.943)    N=32;
01090                 else if(rnd<.97934)  N=34;
01091                 else if(rnd<.99074)  N=33;
01092                 else                 N=36;
01093               }
01094               else            // Cu (29)
01095               {
01096                 if     (rnd<.6917)   N=34;
01097                 else                 N=36;
01098               }
01099             }
01100             else
01101             {
01102               if(i==30)       // Zn
01103               {
01104                 if     (rnd<.486)    N=34;
01105                 else if(rnd<.765)    N=36;
01106                 else if(rnd<.953)    N=38;
01107                 else if(rnd<.994)    N=37;
01108                 else                 N=40;
01109               }
01110               else            // Ga (31)
01111               {
01112                 if     (rnd<.60108)  N=38;
01113                 else                 N=40;
01114               }
01115             }
01116           }
01117         }
01118         else                  // ______ Ge - Mo
01119         {
01120           if     (i<37)       // ...... H - B
01121           {
01122             if   (i<35)
01123             {
01124               if(i==32)       // Ge
01125               {
01126                 if     (rnd<.3594)  N=42;
01127                 else if(rnd<.6360)  N=40;
01128                 else if(rnd<.8484)  N=38;
01129                 else if(rnd<.9256)  N=41;
01130                 else                N=44;
01131               }
01132               else            // Se (34)
01133               {
01134                 if     (rnd>.4961)  N=46;
01135                 else if(rnd<.7378)  N=44;
01136                 else if(rnd<.8274)  N=42;
01137                 else if(rnd<.9148)  N=48;
01138                 else if(rnd<.9911)  N=43;
01139                 else                N=40;
01140               }
01141             }
01142             else
01143             {
01144               if(i==35)       // Br
01145               {
01146                 if     (rnd<.5069)  N=44;
01147                 else                N=46;
01148               }
01149               else            // Kr (36)
01150               {
01151                 if     (rnd<.57)    N=48;
01152                 else if(rnd<.743)   N=50;
01153                 else if(rnd<.859)   N=46;
01154                 else if(rnd<.974)   N=47;
01155                 else if(rnd<.9965)  N=44;
01156                 else                N=42;
01157               }
01158             }
01159           }
01160           else                // ...... Rb - Mo
01161           {
01162             if     (i<40)
01163             {
01164               if(i==37)       // Rb
01165               {
01166                 if     (rnd<.7217)  N=48;
01167                 else                N=50;
01168               }
01169               else            // SR (38)
01170               {
01171                 if     (rnd<.8258)  N=50;
01172                 else if(rnd<.9244)  N=48;
01173                 else if(rnd<.9944)  N=49;
01174                 else                N=46;
01175               }
01176             }
01177             else
01178             {
01179               if(i==40)       // Zr
01180               {
01181                 if     (rnd<.5145)  N=50;
01182                 else if(rnd<.6883)  N=54;
01183                 else if(rnd<.8598)  N=53;
01184                 else if(rnd<.972)   N=51;
01185                 else                N=56;
01186               }
01187               else            // Mo (42)
01188               {
01189                 if     (rnd<.2413)  N=56;
01190                 else if(rnd<.4081)  N=54;
01191                 else if(rnd<.5673)  N=53;
01192                 else if(rnd<.7157)  N=50;
01193                 else if(rnd<.8120)  N=58;
01194                 else if(rnd<.9075)  N=55;
01195                 else                N=52;
01196               }
01197             }
01198           }
01199         }
01200       }
01201     }
01202     else                      // =----= Ru - U
01203     {
01204       if         (i<66)       // ------ Ru - Gd
01205       {
01206         if       (i<54)       // ______ Ru - Te
01207         {
01208           if     (i<49)       // ...... Ru - Cd
01209           {
01210             if   (i<47)
01211             {
01212               if(i==44)       // Ru
01213               {
01214                 if     (rnd<.316)   N=58;
01215                 else if(rnd<.502)   N=60;
01216                 else if(rnd<.673)   N=57;
01217                 else if(rnd<.8)     N=55;
01218                 else if(rnd<.926)   N=56;
01219                 else if(rnd<.9814)  N=52;
01220                 else                N=54;
01221               }
01222               else            // Pd (46)
01223               {
01224                 if     (rnd<.2733)  N=60;
01225                 else if(rnd<.5379)  N=62;
01226                 else if(rnd<.7612)  N=59;
01227                 else if(rnd<.8784)  N=55;
01228                 else if(rnd<.9898)  N=58;
01229                 else                N=56;
01230               }
01231             }
01232             else
01233             {
01234               if(i==47)       // Ag
01235               {
01236                 if(rnd<.51839)      N=60;
01237                 else                N=62;
01238               }
01239               else            // Cd (48)
01240               {
01241                 if     (rnd<.2873)  N=66;
01242                 else if(rnd<.5286)  N=64;
01243                 else if(rnd<.6566)  N=59;
01244                 else if(rnd<.7815)  N=62;
01245                 else if(rnd<.9037)  N=65;
01246                 else if(rnd<.9786)  N=68;
01247                 else if(rnd<.9911)  N=58;
01248                 else                N=60;
01249               }
01250             }
01251           }
01252           else                // ...... In - Te
01253           {
01254             if   (i<51)
01255             {
01256               if(i==49)       // In
01257               {
01258                 if     (rnd<.9577)  N=66;
01259                 else                N=64;
01260               }
01261               else            // Sn (50)
01262               {
01263                 if     (rnd<.3259)  N=70;
01264                 else if(rnd<.5681)  N=68;
01265                 else if(rnd<.7134)  N=66;
01266                 else if(rnd<.7992)  N=69;
01267                 else if(rnd<.8760)  N=67;
01268                 else if(rnd<.9339)  N=74;
01269                 else if(rnd<.9802)  N=72;
01270                 else if(rnd<.9899)  N=62;
01271                 else                N=64;
01272                 //else if(rnd<.9964)  N=64;
01273                 //else                N=65;
01274               }
01275             }
01276             else
01277             {
01278               if(i==51)       // Sb
01279               {
01280                 if     (rnd<.5736)  N=70;
01281                 else                N=72;
01282               }
01283               else            // Te (52)
01284               {
01285                 if     (rnd<.3387)  N=78;
01286                 else if(rnd<.6557)  N=76;
01287                 else if(rnd<.8450)  N=74;
01288                 else if(rnd<.9162)  N=73;
01289                 else if(rnd<.9641)  N=72;
01290                 else if(rnd<.9900)  N=70;
01291                 else if(rnd<.99905) N=71;
01292                 else                N=68;
01293               }
01294             }
01295           }
01296         }
01297         else                // ______ Xe - Gd
01298         {
01299           if     (i<60)     // ...... Xe - B
01300           {
01301             if   (i<57)
01302             {
01303               if(i==54)       // Xe
01304               {
01305                 if     (rnd<.269)   N=78;
01306                 else if(rnd<.533)   N=75;
01307                 else if(rnd<.745)   N=77;
01308                 else if(rnd<.849)   N=80;
01309                 else if(rnd<.938)   N=82;
01310                 else if(rnd<.979)   N=76;
01311                 else if(rnd<.9981)  N=74;
01312                 else if(rnd<.9991)  N=70;
01313                 else                N=72;
01314               }
01315               else            // Ba (56)
01316               {
01317                 if     (rnd<.717)   N=82;
01318                 else if(rnd<.8293)  N=81;
01319                 else if(rnd<.9078)  N=80;
01320                 else if(rnd<.97373) N=79;
01321                 else if(rnd<.99793) N=78;
01322                 else if(rnd<.99899) N=74;
01323                 else                N=76;
01324               }
01325             }
01326             else
01327             {
01328               if(i==57)       // La
01329               {
01330                 if     (rnd<.999098)N=82;
01331                 else                N=81;
01332               }
01333               else            // Ce (58)
01334               {
01335                 if     (rnd<.8843)  N=82;
01336                 else if(rnd<.9956)  N=84;
01337                 else if(rnd<.9981)  N=80;
01338                 else                N=78;
01339               }
01340             }
01341           }
01342           else                // ...... Nd - Gd
01343           {
01344             if   (i<63)
01345             {
01346               if(i==60)       // Nd
01347               {
01348                 if     (rnd<.2713)  N=82;
01349                 else if(rnd<.5093)  N=84;
01350                 else if(rnd<.6812)  N=86;
01351                 else if(rnd<.8030)  N=83;
01352                 else if(rnd<.8860)  N=85;
01353                 else if(rnd<.9436)  N=88;
01354                 else                N=90;
01355               }
01356               else            // Sm (62)
01357               {
01358                 if     (rnd<.267)   N=90;
01359                 else if(rnd<.494)   N=92;
01360                 else if(rnd<.644)   N=85;
01361                 else if(rnd<.782)   N=87;
01362                 else if(rnd<.895)   N=86;
01363                 else if(rnd<.969)   N=88;
01364                 else                N=82;
01365               }
01366             }
01367             else
01368             {
01369               if(i==63)       // Eu
01370               {
01371                 if     (rnd<.522)   N=90;
01372                 else                N=89;
01373               }
01374               else            // Gd (64)
01375               {
01376                 if     (rnd<.2484)  N=94;
01377                 else if(rnd<.4670)  N=96;
01378                 else if(rnd<.6717)  N=92;
01379                 else if(rnd<.8282)  N=93;
01380                 else if(rnd<.9762)  N=91;
01381                 else if(rnd<.9980)  N=90;
01382                 else                N=88;
01383               }
01384             }
01385           }
01386         }
01387       }
01388       else                    // ------ Dy - U
01389       {
01390         if       (i<76)       // ______ Dy - Re
01391         {
01392           if     (i<72)       // ...... Dy - Lu
01393           {
01394             if   (i<70)
01395             {
01396               if(i==66)       // Dy
01397               {
01398                 if     (rnd<.282)   N=98;
01399                 else if(rnd<.537)   N=96;
01400                 else if(rnd<.786)   N=97;
01401                 else if(rnd<.975)   N=95;
01402                 else if(rnd<.9984)  N=94;
01403                 else if(rnd<.9994)  N=92;
01404                 else                N=90;
01405               }
01406               else            // Er (68)
01407               {
01408                 if     (rnd<.3360)  N= 98;
01409                 else if(rnd<.6040)  N=100;
01410                 else if(rnd<.8335)  N= 99;
01411                 else if(rnd<.9825)  N=102;
01412                 else if(rnd<.9986)  N= 96;
01413                 else                N= 94;
01414               }
01415             }
01416             else
01417             {
01418               if(i==70)       // Yb
01419               {
01420                 if     (rnd<.3180)  N=104;
01421                 else if(rnd<.5370)  N=102;
01422                 else if(rnd<.6982)  N=103;
01423                 else if(rnd<.8412)  N=101;
01424                 else if(rnd<.9682)  N=106;
01425                 else if(rnd<.9987)  N=100;
01426                 else                N= 98;
01427               }
01428               else            // Lu (71)
01429               {
01430                 if     (rnd<.9741)  N=104;
01431                 else                N=105;
01432               }
01433             }
01434           }
01435           else                // ...... Hf - Re
01436           {
01437             if   (i<74)
01438             {
01439               if(i==72)       // Hf
01440               {
01441                 if     (rnd<.35100) N=108;
01442                 else if(rnd<.62397) N=106;
01443                 else if(rnd<.81003) N=105;
01444                 else if(rnd<.94632) N=107;
01445                 else if(rnd<.99838) N=104;
01446                 else                N=102;
01447               }
01448               else            // Ta (73)
01449               {
01450                 if(rnd<.99988) N=108;
01451                 else           N=107;
01452               }
01453             }
01454             else
01455             {
01456               if(i==74)       // W
01457               {
01458                 if     (rnd<.307)   N=110;
01459                 else if(rnd<.593)   N=112;
01460                 else if(rnd<.856)   N=108;
01461                 else if(rnd<.9988)  N=109;
01462                 else                N=106;
01463               }
01464               else            // Re (75)
01465               {
01466                 if     (rnd<.626)   N=112;
01467                 else                N=110;
01468               }
01469             }
01470           }
01471         }
01472         else                  // ______ Os - U
01473         {
01474           if     (i<81)       // ...... Os - Hg
01475           {
01476             if     (i<78)
01477             {
01478               if(i==76)       // Os
01479               {
01480                 if     (rnd<.410)   N=116;
01481                 else if(rnd<.674)   N=114;
01482                 else if(rnd<.835)   N=113;
01483                 else if(rnd<.968)   N=112;
01484                 else if(rnd<.984)   N=111;
01485                 else if(rnd<.9998)  N=110;
01486                 else                N=108;
01487               }
01488               else            // Ir (77)
01489               {
01490                 if     (rnd<.627)   N=116;
01491                 else                N=114;
01492               }
01493             }
01494             else
01495             {
01496               if(i==78)       // Pt
01497               {
01498                 if     (rnd<.338)   N=117;
01499                 else if(rnd<.667)   N=116;
01500                 else if(rnd<.920)   N=118;
01501                 else if(rnd<.992)   N=120;
01502                 else if(rnd<.9999)  N=114;
01503                 else                N=112;
01504               }
01505               else            // Hg (80)
01506               {
01507                 if     (rnd<.2986)  N=122;
01508                 else if(rnd<.5296)  N=120;
01509                 else if(rnd<.6983)  N=119;
01510                 else if(rnd<.8301)  N=121;
01511                 else if(rnd<.9298)  N=118;
01512                 else if(rnd<.9985)  N=124;
01513                 else                N=116;
01514               }
01515             }
01516           }
01517           else                // ...... Tl - U
01518           {
01519             if        (i<92)
01520             {
01521               if     (i==81)  // Tl
01522               {
01523                 if     (rnd<.70476) N=124;
01524                 else                N=122;
01525               }
01526               else            // Pb (82)
01527               {
01528                 if     (rnd<.524)   N=126;
01529                 else if(rnd<.765)   N=124;
01530                 else if(rnd<.986)   N=125;
01531                 else                N=122;
01532               }
01533             }
01534             else              // U (92)
01535             {
01536                 if     (rnd<.992745)N=146;
01537                 else if(rnd<.999945)N=143;
01538                 else                N=142;
01539             }
01540           }
01541         }
01542       }
01543     }
01544   }
01545   else if(N<0)
01546   {
01547     N=-N;
01548     G4cerr<<"---Worn---G4QIsotope::RandomizeNeutrons:UnstableElement's used Z="<<i<<G4endl;
01549   }
01550 #ifdef debug
01551   G4cout<<"G4QIsotope::RandomizeNeutrons: Z="<<i<<", N="<<N<<G4endl;
01552 #endif
01553   return N;
01554 }
01555 
01556 // Returns the input index (if it is >0 & unique) or theFirstFreeIndex (<=0 or nonunique)
01557 G4int G4QIsotope::InitElement(G4int Z, G4int index, // Ret: -1 - Empty, -2 - Wrong (sum>1)
01558                               vector<pair<G4int,G4double>*>* abund)
01559 {
01560   G4int I=abund->size();
01561 #ifdef debug
01562   G4cout<<"G4QIsotope::InitElement: called with I="<<I<<" pairs,Z="<<Z<<",i="<<indexG4endl;
01563 #endif
01564   if(I<=0)
01565   {
01566     G4cerr<<"--Worning--G4QIsotope::InitEl:(-1)0VectorOfNewEl,Z="<<Z<<",i="<<index<<G4endl;
01567     return -2;
01568   }
01569   if(IsDefined(Z,index))              // This index is already defined
01570   {
01571     G4cerr<<"-Worning-G4QIsotope::InitEl:VONewEl,Z="<<Z<<",ind="<<index<<" exists"<<G4endl;
01572     return index;
01573   }
01574   G4int ZInd=1000*index+Z;            // Fake Z increased by the UserDefinedIndex
01575   vector<pair<G4int,G4double>*>*A = new vector<pair<G4int,G4double>*>;
01576   vector<pair<G4int,G4double>*>*S = new vector<pair<G4int,G4double>*>;
01577   vector<pair<G4int,G4double>*>*C = new vector<pair<G4int,G4double>*>;
01578 #ifdef debug
01579   G4cout<<"G4QIsotope::InitElement: A & S & C vectors are alocated"<<G4endl;
01580 #endif
01581   G4double sumAbu=0;                  // Summ of abbundancies
01582   for(G4int j=0; j<I; j++)
01583   {
01584     G4int N=(*abund)[j]->first;
01585     G4double abu=(*abund)[j]->second;
01586 #ifdef debug
01587     G4cout<<"G4QIsotope::InitElement: pair#"<<j<<", N="<<N<<", abund="<<abu<<G4endl;
01588 #endif
01589     sumAbu+=abu;
01590     if(j==I-1.)
01591     {
01592       if(fabs(sumAbu-1.)>.00001)
01593       {
01594         G4cerr<<"--Worning--G4QIsotope::InitEl:maxSum="<<sumAbu<<" is fixed to 1."<<G4endl;
01595         abu+=1.-sumAbu;
01596         sumAbu=1.;
01597       }
01598       else if(sumAbu-abu>1.)
01599       {
01600         G4cerr<<"--Worning--G4QIsotope::InitEl:(-2)WrongAbund,Z="<<Z<<",i="<<index<<G4endl;
01601         for(G4int k=0; k<I-1; k++)
01602         {
01603           delete (*A)[k];
01604           delete (*S)[k];
01605           delete (*C)[k];
01606         }
01607         delete A; 
01608         delete S;
01609         delete C;
01610         return -2;
01611       }
01612 #ifdef debug
01613     G4cout<<"G4QIsotope::InitElement:TheLastIsChecked it isOK or coredTo "<<sumAbu<<G4endl;
01614 #endif
01615     }
01616     pair<G4int,G4double>* abP= new pair<G4int,G4double>(N,abu);
01617     A->push_back(abP); // @@ Valgrind thinks that it is not deleted (?) (Line 703)
01618     pair<G4int,G4double>* saP= new pair<G4int,G4double>(N,sumAbu);
01619     S->push_back(saP); // @@ Valgrind thinks that it is not deleted (?) (Line 713)
01620     pair<G4int,G4double>* csP= new pair<G4int,G4double>(N,0.);
01621     C->push_back(csP); // @@ Valgrind thinks that it is not deleted (?) (Line 723)
01622 #ifdef debug
01623     G4cout<<"G4QIsotope::InitElement: A & S & C are filled nP="<<C->size()<<G4endl;
01624 #endif
01625   }
01626   pair<G4int,vector<pair<G4int,G4double>*>*>* newAP=
01627                                     new pair<G4int,vector<pair<G4int,G4double>*>*>(ZInd,A);
01628   newElems.push_back(newAP);
01629   pair<G4int,vector<pair<G4int,G4double>*>*>* newSA=
01630                                     new pair<G4int,vector<pair<G4int,G4double>*>*>(ZInd,S);
01631   newSumAb.push_back(newSA);
01632   pair<G4int,vector<pair<G4int,G4double>*>*>* newCP=
01633                                     new pair<G4int,vector<pair<G4int,G4double>*>*>(ZInd,C);
01634   newIsoCS.push_back(newCP);
01635 #ifdef debug
01636   G4cout<<"G4QIsotope::InitElement: newElems & newSumAb & newIsoCS are filled "<<G4endl;
01637 #endif
01638   return index;
01639 }
01640 
01641 // The highest index defined for Element with Z (Index>0 correspondToUserDefinedElements)
01642 G4int G4QIsotope::GetLastIndex(G4int Z) // Returns theLastDefinedIndex (if onlyNatural: =0)
01643 {
01644 #ifdef debug
01645   G4cout<<"G4QIsotope::GetLastIndex is called Z="<<Z<<G4endl;
01646 #endif
01647   G4int mind=0;                       // Prototype of the maximum existing index for this Z
01648   G4int nE=newElems.size();           // A number of definitions in the newElements vector
01649   if(nE) for(G4int i=0; i<nE; i++)    // LOOP over new UserDefinedElements
01650   {
01651     G4int zin=newElems[i]->first;
01652     G4int zi=zin%1000;                // Existing Z
01653     G4int in=zin/1000;                // Existing index
01654     if(Z==zi && in>mind) mind=in;     // maximum index for this Z
01655   }
01656   return mind;
01657 }
01658 
01659 // Indices can have differen numbers (not 1,2,3,...) & in different sequences (9,3,7,...)
01660 G4bool G4QIsotope::IsDefined(G4int Z, G4int Ind) // Ind is an index to be found (true)
01661 {
01662 #ifdef debug
01663   G4cout<<"G4QIsotope::IsDefined is called Z="<<Z<<", I="<<Ind<<G4endl;
01664 #endif
01665   if(Ind<=0)
01666   {
01667     if(Ind<0) G4cerr<<"-W-G4QIsotope::IsDefined: Z="<<Z<<", Ind="<<Ind<<" < 0->=0"<<G4endl;
01668     return true;                      // to avoid definition with the negative index
01669   }
01670   G4int nE=newElems.size();           // A number of definitions in the newElements vector
01671   if(nE) for(G4int i=0; i<nE; i++)    // LOOP over new UserDefinedElements
01672   {
01673     G4int zin=newElems[i]->first;
01674     G4int zi=zin%1000;                // Existing Z
01675     G4int in=zin/1000;                // Existing index
01676     if(Z==zi && Ind==in) return true;  // The index for the element Z is found
01677   }
01678   return false;                       // The index for the element Z is not found
01679 }
01680 
01681 // A#ofNeutrons in theElement with Z & UserDefIndex. Universal for Nat(index=0) & UserDefEl
01682 G4int G4QIsotope::GetNeutrons(G4int Z, G4int index) // If theElem doesn't exist, returns <0
01683 {
01684 #ifdef debug
01685   G4cout<<"G4QIsotope::GetNeutrons is called Z="<<Z<<", index="<<index<<G4endl;
01686 #endif
01687   // To reduce the code, but make the member function a bit slower, one can use for natural
01688   // isotopes the same algorithm as for the newElements, splitting the natElements Vector
01689   if(!index) return RandomizeNeutrons(Z); // @@ Fast decision for the natural isotopes
01690   else if(index<0)
01691   {
01692     G4cerr<<"---Worning---G4QIsotope::GetNeutrons:(-2) Negative Index i="<<index<<G4endl;
01693     return -2;
01694   }
01695   // For the positive index tries to randomize the newUserDefinedElement
01696   G4bool found=false;                 // Prototype of the"ZWithTheSameIndex is found" event
01697   G4int nE=newElems.size();           // A number of definitions in the newElements Vector
01698   G4int i=0;
01699   if(nE) for(i=0; i<nE; i++)
01700   {
01701     G4int zin=newElems[i]->first;
01702     G4int in=zin/1000;                // Existing index
01703     G4int zi=zin%1000;                // Existing Z
01704     if(Z==zi && in==index)
01705     {
01706       found=true;                     // The newElement with the same Z & index is found
01707       break;                          // Finish the search and quit the loop
01708     }
01709   }
01710   if(!found)
01711   {
01712     G4cerr<<"--Worning--G4QIsotope::GetNeutrons:(-1) NotFound Z="<<Z<<",i="<<index<<G4endl;
01713     return -1;
01714   }
01715   vector<pair<G4int,G4double>*>* abu = newSumAb[i]->second;
01716   G4int nn = abu->size();             // A#Of UserDefinedIsotopes for the newElement
01717   if(nn>0)
01718   {
01719     if(nn==1) return (*abu)[0]->first;
01720     else
01721     {
01722       G4double rnd=G4UniformRand();
01723       G4int j=0;
01724       for(j=0; j<nn; j++) if ( rnd < (*abu)[j]->second ) break;
01725       if(j>=nn) j=nn-1;
01726       return (*abu)[j]->first;
01727     }
01728   }
01729   else
01730   {
01731     G4cerr<<"--Worning--G4QIsotope::GetNeutrons:(-3) Empty Z="<<Z<<",i="<<index<<G4endl;
01732     return -3;
01733   }
01734 }
01735 
01736 // Get a pointer to the vector of pairs(N,CrosSec), where N is used to calculate CrosSec
01737 vector<pair<G4int,G4double>*>* G4QIsotope::GetCSVector(G4int Z, G4int index)
01738 {
01739 #ifdef debug
01740   G4cout<<"G4QIsotope::GetCSVector is called"<<G4endl;
01741 #endif
01742   if(index<0)
01743   {
01744     G4cerr<<"---Worning---G4QIsotope::GetSCVector:(-1) Negative Index i="<<index<<G4endl;
01745     return 0;
01746   }
01747   else if(!index) return natIsoCrosS[Z];
01748   // For the positive index tries to find the newUserDefinedElement
01749   G4bool found=false;                 // Prototype of the"ZWithTheSameIndex is found" event
01750   G4int nE=newIsoCS.size();           // A number of definitions in the newElements Vector
01751   G4int i=0;
01752   if(nE) for(i=0; i<nE; i++)
01753   {
01754     G4int zin=newIsoCS[i]->first;
01755     G4int in=zin/1000;                // Existing index
01756     G4int zi=zin%1000;                // Existing Z
01757     if(Z==zi && in==index)
01758     {
01759       found=true;                     // The newElement with the same Z & index is found
01760       break;                          // Finish the search and quit the loop
01761     }
01762   }
01763   if(!found)
01764   {
01765     G4cerr<<"--Worning--G4QIsotope::GetSCVector:(-2) NotFound Z="<<Z<<",i="<<index<<G4endl;
01766     return 0;
01767   }
01768   return newIsoCS[i]->second;
01769 }
01770 
01771 // Get a pointer to the vector of pairs(N,IntAbundancy) for the element with Z
01772 vector<pair<G4int,G4double>*>* G4QIsotope::GetAbuVector(G4int Z, G4int index)
01773 {
01774 #ifdef debug
01775   G4cout<<"G4QIsotope::GetAbuVector is called"<<G4endl;
01776 #endif
01777   if(index<0)
01778   {
01779     G4cerr<<"---Worning---G4QIsotope::GetAbuVector:(-1) Negative Index i="<<index<<G4endl;
01780     return 0;
01781   }
01782   else if(!index) return natElements[Z];
01783   // For the positive index tries to find the newUserDefinedElement
01784   G4bool found=false;                 // Prototype of the"ZWithTheSameIndex is found" event
01785   G4int nE=newElems.size();           // A number of definitions in the newElements Vector
01786   G4int i=0;
01787   if(nE) for(i=0; i<nE; i++)
01788   {
01789     G4int zin=newElems[i]->first;
01790     G4int in=zin/1000;                // Existing index
01791     G4int zi=zin%1000;                // Existing Z
01792     if(Z==zi && in==index)
01793     {
01794       found=true;                     // The newElement with the same Z & index is found
01795       break;                          // Finish the search and quit the loop
01796     }
01797   }
01798   if(!found)
01799   {
01800     G4cerr<<"--Worning--G4QIsotope::GetAbuVector:(-2)NotFound Z="<<Z<<",i="<<index<<G4endl;
01801     return 0;
01802   }
01803   return newElems[i]->second;
01804 }
01805 
01806 // Get a pointer to the vector of pairs(N,SumAbundancy) for the element with Z
01807 vector<pair<G4int,G4double>*>* G4QIsotope::GetSumAVector(G4int Z, G4int index)
01808 {
01809 #ifdef debug
01810   G4cout<<"G4QIsotope::GetSumAVector is called"<<G4endl;
01811 #endif
01812   if(index<0)
01813   {
01814     G4cerr<<"---Worning---G4QIsotope::GetSumAVector:(-1) Negative Index i="<<index<<G4endl;
01815     return 0;
01816   }
01817   else if(!index) return natSumAbund[Z];
01818   // For the positive index tries to find the newUserDefinedElement
01819   G4bool found=false;                 // Prototype of the"ZWithTheSameIndex is found" event
01820   G4int nE=newSumAb.size();           // A number of definitions in the newElements Vector
01821   G4int i=0;
01822   if(nE) for(i=0; i<nE; i++)
01823   {
01824     G4int zin=newSumAb[i]->first;
01825     G4int in=zin/1000;                // Existing index
01826     G4int zi=zin%1000;                // Existing Z
01827     if(Z==zi && in==index)
01828     {
01829       found=true;                     // The newElement with the same Z & index is found
01830       break;                          // Finish the search and quit the loop
01831     }
01832   }
01833   if(!found)
01834   {
01835     G4cerr<<"-Worning-G4QIsotope::GetSumAVector:(-2)Not Found Z="<<Z<<",i="<<index<<G4endl;
01836     return 0;
01837   }
01838   return newSumAb[i]->second;
01839 }
01840 
01841 // Calculates the mean Cross Section for the initialized Element(ind=0 Nat,ind>0 UserDef)
01842 G4double G4QIsotope::GetMeanCrossSection(G4int Z, G4int index)
01843 {
01844   vector<pair<G4int,G4double>*>* ab;
01845   vector<pair<G4int,G4double>*>* cs;
01846 #ifdef ppdebug
01847   G4cout<<"G4QIsotope::GetMeanCrossSection is called"<<G4endl;
01848 #endif
01849   if(index<0)
01850   {
01851     G4cerr<<"---Worning---G4QIsotope::GetMeanCS:(-1) Negative Index i="<<index<<G4endl;
01852     return -1.;
01853   }
01854   else if(!index)           // =-------=> Natural Abundancies for Isotopes of the Element
01855   {
01856 #ifdef ppdebug
01857     G4cout<<"G4QIsotope::GetMeanCrossSection: Nat Abundance, Z="<<Z<<G4endl;
01858 #endif
01859     ab=natElements[Z];
01860     cs=natIsoCrosS[Z];
01861   }
01862   else                      // =------=> UserDefinedAbundancies for Isotopes of theElement
01863   {
01864 #ifdef ppdebug
01865     G4cout<<"G4QIsotope::GetMeanCrossSection: Art Abund, Z="<<Z<<",ind="<<index<<G4endl;
01866 #endif
01867     // For the positive index tries to find the newUserDefinedElement
01868     G4bool found=false;               // Prototype of the"ZWithTheSameIndex is found" event
01869     G4int nE=newIsoCS.size();         // A number of definitions in the newElements Vector
01870     G4int i=0;
01871     if(nE) for(i=0; i<nE; i++)
01872     {
01873       G4int zin=newIsoCS[i]->first;
01874       G4int in=zin/1000;              // Existing index
01875       G4int zi=zin%1000;              // Existing Z
01876       if(Z==zi && in==index)
01877       {
01878         found=true;                   // The newElement with the same Z & index is found
01879         break;                        // Finish the search and quit the loop
01880       }
01881     }
01882     if(!found)
01883     {
01884       G4cerr<<"--Worning--G4QIsotope::GetMeanCS:(-2) NotFound Z="<<Z<<",i="<<index<<G4endl;
01885       return -2.;
01886     }
01887     ab=newElems[i]->second;
01888     cs=newIsoCS[i]->second;
01889   }
01890   G4int nis=ab->size();
01891   //G4double last=0.;
01892   if(!nis)
01893   {
01894     G4cerr<<"--Worning--G4QIsotope::GetMeanCS:(-3) Empty Z="<<Z<<",i="<<index<<G4endl;
01895     return -3.;
01896   }
01897   else
01898   {
01899     G4double sum=0.;
01900     for(G4int j=0; j<nis; j++)
01901     {
01902       G4double cur=(*ab)[j]->second;
01903       //G4double abunda=cur-last;
01904       //last=cur;
01905 #ifdef ppdebug
01906       G4cout<<"G4QIsot::GetMeanCS:j="<<j<<",ab="<<cur<<",CS="<<(*cs)[j]->second<<G4endl;
01907 #endif
01908       //sum+=abunda * (*cs)[j]->second;
01909       sum+=cur * (*cs)[j]->second;
01910     }
01911     return sum;
01912   }
01913 }
01914 
01915 // Randomize A#OfNeutrons in the Isotope weighted by theAbubdancies and theCrossSections
01916 G4int G4QIsotope::GetCSNeutrons(G4int Z, G4int index)
01917 {
01918   vector<pair<G4int,G4double>*>* ab;
01919   vector<pair<G4int,G4double>*>* cs;
01920 #ifdef debug
01921   G4cout<<"G4QIsotope::GetCSNeutrons is called"<<G4endl;
01922 #endif
01923   if(index<0)
01924   {
01925     G4cerr<<"---Worning---G4QIsotope::GetCSNeutrons:(-1) Negative Index i="<<index<<G4endl;
01926     return -1;
01927   }
01928   else if(!index)           // =---------=> Natural Abundancies for Isotopes of the Element
01929   {
01930     ab=natElements[Z];
01931     cs=natIsoCrosS[Z];
01932   }
01933   else                      // =-------=> UserDefinedAbundancies for Isotopes of theElement
01934   {
01935     // For the positive index tries to find the newUserDefinedElement
01936     G4bool found=false;               // Prototype of the"ZWithTheSameIndex is found" event
01937     G4int nE=newIsoCS.size();         // A number of definitions in the newElements Vector
01938     G4int i=0;
01939     if(nE) for(i=0; i<nE; i++)
01940     {
01941       G4int zin=newIsoCS[i]->first;
01942       G4int in=zin/1000;              // Existing index
01943       G4int zi=zin%1000;              // Existing Z
01944       if(Z==zi && in==index)
01945       {
01946         found=true;                   // The newElement with the same Z & index is found
01947         break;                        // Finish the search and quit the loop
01948       }
01949     }
01950     if(!found)
01951     {
01952       G4cerr<<"--Worning--G4QIsotope::GetCSNeut:(-2) NotFound Z="<<Z<<",i="<<index<<G4endl;
01953       return -2;
01954     }
01955     ab=newElems[i]->second;
01956     cs=newIsoCS[i]->second;
01957   }
01958   G4int nis=ab->size();
01959   G4double last=0.;
01960   if(!nis)
01961   {
01962     G4cerr<<"--Worning--G4QIsotope::GetCSNeutrons:(-3) Empty Z="<<Z<<",i="<<index<<G4endl;
01963     return -3;
01964   }
01965   else
01966   {
01967     G4double sum=0.;
01968     vector<G4double> scs(nis);
01969     for(G4int j=0; j<nis; j++)
01970     {
01971       G4double cur=(*ab)[j]->second;
01972       G4double abunda=cur-last;
01973       last=cur;
01974       sum+=abunda * (*cs)[j]->second;;
01975       scs.push_back(sum);
01976     }
01977     G4double rnd=sum*G4UniformRand();
01978     sum=0;
01979     G4int k=0;
01980     if(nis>1) for(k=0; k<nis; k++) if(rnd<scs[k]) break;
01981     return (*ab)[k]->first;
01982   }
01983 }

Generated on Mon May 27 17:49:36 2013 for Geant4 by  doxygen 1.4.7