G4RadioactiveDecay.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 // MODULE:              G4RadioactiveDecay.cc
00027 //
00028 // Author:              F Lei & P R Truscott
00029 // Organisation:        DERA UK
00030 // Customer:            ESA/ESTEC, NOORDWIJK
00031 // Contract:            12115/96/JG/NL Work Order No. 3
00032 //
00033 // Documentation avaialable at http://www.space.dera.gov.uk/space_env/rdm.html
00034 //   These include:
00035 //       User Requirement Document (URD)
00036 //       Software Specification Documents (SSD)
00037 //       Software User Manual (SUM)
00038 //       Technical Note (TN) on the physics and algorithms
00039 //
00040 //    The test and example programs are not included in the public release of 
00041 //    G4 but they can be downloaded from
00042 //      http://www.space.qinetiq.com/space_env/rdm.html
00043 // 
00044 // CHANGE HISTORY
00045 // --------------
00046 // 03 Oct  2012, V. Ivanchenko removed internal table for mean free path 
00047 //                             similar to what is done for as G4Decay
00048 // 10 July 2012, L. Desorgher
00049 //                      -In LoadDecayTable:  Add LoadedNuclei.push_back(theParentNucleus.GetParticleName());
00050 //                      also for the case where user data files are used. Correction for bug
00051 //                      1324. Changes proposed by Joa L.
00052 //
00053 //
00054 // 01 May 2012, L. Desorgher
00055 //                      -Force the reading of user file to theIsotopeTable
00056 //                      -Merge the development by Fan Lei for activation computation
00057 //
00058 // 17 Oct 2011, L. Desorgher
00059 //                      -Add possibility for the user to load its own decay file.
00060 //                      -Set halflifethreshold negative by default to allow the tracking of all
00061 //                         excited nuclei resulting from a radioactive decay
00062 //
00063 // 01 June 2011, M. Kelsey -- Add directional biasing interface to allow for
00064 //              "collimation" of decay daughters.
00065 // 16 February 2006, V.Ivanchenko fix problem in IsApplicable connected with
00066 //            8.0 particle design
00067 // 18 October 2002, F. Lei
00068 //            in the case of beta decay, added a check of the end-energy 
00069 //            to ensure it is > 0.
00070 //            ENSDF occationally have beta decay entries with zero energies
00071 //
00072 // 27 Sepetember 2001, F. Lei
00073 //            verboselevel(0) used in constructor
00074 //
00075 // 01 November 2000, F.Lei
00076 //            added " ee = e0 +1. ;" as line 763
00077 //            tagged as "radiative_decay-V02-00-02"              
00078 // 28 October 2000, F Lei 
00079 //            added fast beta decay mode. Many files have been changed.
00080 //            tagged as "radiative_decay-V02-00-01"
00081 //
00082 // 25 October 2000, F Lei, DERA UK
00083 //            1) line 1185 added 'const' to work with tag "Track-V02-00-00"
00084 //            tagged as "radiative_decay-V02-00-00"
00085 // 14 April 2000, F Lei, DERA UK
00086 // 0.b.4 release. Changes are:
00087 //            1) Use PhotonEvaporation instead of DiscreteGammaDeexcitation
00088 //            2) VR: Significant efficiency improvement
00089 // 
00090 // 29 February 2000, P R Truscott, DERA UK
00091 // 0.b.3 release.
00092 //
00094 //
00095 #include "G4RadioactiveDecay.hh"
00096 #include "G4RadioactiveDecaymessenger.hh"
00097 
00098 #include "G4SystemOfUnits.hh"
00099 #include "G4DynamicParticle.hh"
00100 #include "G4DecayProducts.hh"
00101 #include "G4DecayTable.hh"
00102 #include "G4PhysicsLogVector.hh"
00103 #include "G4ParticleChangeForRadDecay.hh"
00104 #include "G4ITDecayChannel.hh"
00105 #include "G4BetaMinusDecayChannel.hh"
00106 #include "G4BetaPlusDecayChannel.hh"
00107 #include "G4KshellECDecayChannel.hh"
00108 #include "G4LshellECDecayChannel.hh"
00109 #include "G4MshellECDecayChannel.hh"
00110 #include "G4AlphaDecayChannel.hh"
00111 #include "G4VDecayChannel.hh"
00112 #include "G4RadioactiveDecayMode.hh"
00113 #include "G4Ions.hh"
00114 #include "G4IonTable.hh"
00115 #include "G4RIsotopeTable.hh"
00116 #include "G4BetaDecayType.hh"
00117 #include "G4BetaDecayCorrections.hh"
00118 #include "Randomize.hh"
00119 #include "G4LogicalVolumeStore.hh"
00120 #include "G4NuclearLevelManager.hh"
00121 #include "G4NuclearLevelStore.hh"
00122 #include "G4ThreeVector.hh"
00123 #include "G4Electron.hh"
00124 #include "G4Positron.hh"
00125 #include "G4Neutron.hh"
00126 #include "G4Gamma.hh"
00127 #include "G4Alpha.hh"
00128 
00129 #include "G4HadTmpUtil.hh"
00130 #include "G4HadronicProcessType.hh"
00131 #include "G4LossTableManager.hh"
00132 #include "G4VAtomDeexcitation.hh"
00133 
00134 #include <vector>
00135 #include <sstream>
00136 #include <algorithm>
00137 #include <fstream>
00138 
00139 using namespace CLHEP;
00140 
00141 const G4double G4RadioactiveDecay::levelTolerance =2.0*keV;
00142 const G4ThreeVector G4RadioactiveDecay::origin(0.,0.,0.);
00143 
00144 G4RadioactiveDecay::G4RadioactiveDecay(const G4String& processName)
00145  : G4VRestDiscreteProcess(processName, fDecay), HighestValue(20.0),
00146    isInitialised(false), forceDecayDirection(0.,0.,0.), 
00147    forceDecayHalfAngle(0.*deg), verboseLevel(0)
00148 {
00149 #ifdef G4VERBOSE
00150   if (GetVerboseLevel()>1) {
00151     G4cout <<"G4RadioactiveDecay constructor    Name: ";
00152     G4cout <<processName << G4endl;   }
00153 #endif
00154 
00155   SetProcessSubType(fRadioactiveDecay);
00156 
00157   theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
00158   theIsotopeTable              = new G4RIsotopeTable();
00159   pParticleChange              = &fParticleChangeForRadDecay;
00160 
00161   // Now register the Isotope table with G4IonTable.
00162   G4IonTable *theIonTable =
00163     (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
00164   G4VIsotopeTable *aVirtualTable = theIsotopeTable;
00165   theIonTable->RegisterIsotopeTable(aVirtualTable);
00166 
00167   // Reset the contents of the list of nuclei for which decay scheme data
00168   // have been loaded.
00169   LoadedNuclei.clear();
00170 
00171   //
00172   //Reset the list of user define data file
00173   //
00174   theUserRadioactiveDataFiles.clear();
00175 
00176   //
00177   //
00178   // Apply default values.
00179   //
00180   NSourceBin  = 1;
00181   SBin[0]     = 0.* s;
00182   SBin[1]     = 1.* s;
00183   SProfile[0] = 1.;
00184   SProfile[1] = 0.;
00185   NDecayBin   = 1;
00186   DBin[0]     = 0. * s ;
00187   DBin[1]     = 1. * s;
00188   DProfile[0] = 1.;
00189   DProfile[1] = 0.;
00190   decayWindows[0] = 0;
00191   G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
00192   theRadioactivityTables.push_back(rTable);
00193   NSplit      = 1;
00194   AnalogueMC  = true ;
00195   FBeta       = false ;
00196   BRBias      = true ;
00197   applyICM    = true ;
00198   applyARM    = true ;
00199   halflifethreshold = -1.*second;
00200   //
00201   // RDM applies to xall logical volumes as default
00202   isAllVolumesMode=true;
00203   SelectAllVolumes();
00204 }
00205 
00206 
00207 G4RadioactiveDecay::~G4RadioactiveDecay()
00208 {
00209   delete theRadioactiveDecaymessenger;
00210 }
00211 
00212 
00213 G4bool
00214 G4RadioactiveDecay::IsApplicable(const G4ParticleDefinition& aParticle)
00215 {
00216   // All particles, other than G4Ions, are rejected by default
00217   if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
00218   if (aParticle.GetParticleName() == "GenericIon") {
00219     return true;
00220   } else if (!(aParticle.GetParticleType() == "nucleus")
00221              || aParticle.GetPDGLifeTime() < 0. ) {
00222     return false;
00223   }
00224 
00225   // Determine whether the nuclide falls into the correct A and Z range
00226 
00227   G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
00228   G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
00229   if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
00230     {return false;}
00231   else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
00232     {return false;}
00233   return true;
00234 }
00235 
00236 
00237 G4bool G4RadioactiveDecay::IsLoaded(const G4ParticleDefinition &aParticle)
00238 {
00239   // Check whether the radioactive decay data on the ion have already been
00240   // loaded
00241 
00242   return std::binary_search(LoadedNuclei.begin(),
00243                             LoadedNuclei.end(),
00244                             aParticle.GetParticleName());
00245 }
00246 
00247 
00248 void G4RadioactiveDecay::SelectAVolume(const G4String aVolume)
00249 {
00250   G4LogicalVolumeStore* theLogicalVolumes;
00251   G4LogicalVolume *volume;
00252   theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
00253   for (size_t i = 0; i < theLogicalVolumes->size(); i++){
00254     volume=(*theLogicalVolumes)[i];
00255     if (volume->GetName() == aVolume) {
00256       ValidVolumes.push_back(aVolume);
00257       std::sort(ValidVolumes.begin(), ValidVolumes.end());
00258       // sort need for performing binary_search
00259 #ifdef G4VERBOSE
00260       if (GetVerboseLevel()>0)
00261         G4cout << " RDM Applies to : " << aVolume << G4endl; 
00262 #endif
00263     }else if(i ==  theLogicalVolumes->size())
00264       {
00265         G4cerr << "SelectAVolume: "<<aVolume << " is not a valid logical volume name"<< G4endl; 
00266       }
00267   }
00268 }
00269 
00270 
00271 void G4RadioactiveDecay::DeselectAVolume(const G4String aVolume)
00272 {
00273   G4LogicalVolumeStore *theLogicalVolumes;
00274   G4LogicalVolume *volume;
00275   theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
00276   for (size_t i = 0; i < theLogicalVolumes->size(); i++){
00277     volume=(*theLogicalVolumes)[i];
00278     if (volume->GetName() == aVolume) {
00279       std::vector<G4String>::iterator location;
00280       location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
00281       if (location != ValidVolumes.end()) {
00282         ValidVolumes.erase(location);
00283         std::sort(ValidVolumes.begin(), ValidVolumes.end());
00284         isAllVolumesMode =false;
00285       } else {
00286         G4cerr << " DeselectVolume:" << aVolume << " is not in the list"<< G4endl; 
00287       }   
00288 #ifdef G4VERBOSE
00289       if (GetVerboseLevel()>0)
00290         G4cout << " DeselectVolume: " << aVolume << " is removed from list"<<G4endl; 
00291 #endif
00292     } else if (i ==  theLogicalVolumes->size()) {
00293       G4cerr << " DeselectVolume:" << aVolume
00294              << "is not a valid logical volume name"<< G4endl; 
00295     }
00296   }
00297 }
00298 
00299 
00300 void G4RadioactiveDecay::SelectAllVolumes() 
00301 {
00302   G4LogicalVolumeStore *theLogicalVolumes;
00303   G4LogicalVolume *volume;
00304   theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
00305   ValidVolumes.clear();
00306 #ifdef G4VERBOSE
00307   if (GetVerboseLevel()>0)
00308     G4cout << " RDM Applies to all Volumes"  << G4endl;
00309 #endif
00310   for (size_t i = 0; i < theLogicalVolumes->size(); i++){
00311     volume=(*theLogicalVolumes)[i];
00312     ValidVolumes.push_back(volume->GetName());    
00313 #ifdef G4VERBOSE
00314     if (GetVerboseLevel()>0)
00315       G4cout << "         RDM Applies to Volume "  << volume->GetName() << G4endl;
00316 #endif
00317   }
00318   std::sort(ValidVolumes.begin(), ValidVolumes.end());
00319   // sort needed in order to allow binary_search
00320   isAllVolumesMode=true;
00321 }
00322 
00323 
00324 void G4RadioactiveDecay::DeselectAllVolumes() 
00325 {
00326   ValidVolumes.clear();
00327   isAllVolumesMode=false;
00328 #ifdef G4VERBOSE
00329   if (GetVerboseLevel()>0)
00330     G4cout << " RDM removed from all volumes" << G4endl; 
00331 #endif
00332 }
00333 
00334 
00335 G4bool
00336 G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition& aParticle)
00337 {
00338   // Check whether the radioactive decay rates table for the ion has already
00339   // been calculated.
00340   G4String aParticleName = aParticle.GetParticleName();
00341   for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
00342     if (theDecayRateTableVector[i].GetIonName() == aParticleName)
00343       return true;
00344   }
00345   return false;
00346 }
00347 
00348 // GetDecayRateTable
00349 // retrieve the decayratetable for the specified aParticle
00350 
00351 void
00352 G4RadioactiveDecay::GetDecayRateTable(const G4ParticleDefinition& aParticle)
00353 {
00354   G4String aParticleName = aParticle.GetParticleName();
00355 
00356   for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
00357     if (theDecayRateTableVector[i].GetIonName() == aParticleName) {
00358       theDecayRateVector = theDecayRateTableVector[i].GetItsRates();
00359     }
00360   }
00361 #ifdef G4VERBOSE
00362   if (GetVerboseLevel() > 0) {
00363     G4cout << "The DecayRate Table for "
00364            << aParticleName << " is selected." <<  G4endl;
00365   }
00366 #endif
00367 }
00368 
00369 // GetTaoTime performs the convolution of the source time profile function
00370 // with the decay constants in the decay chain. 
00371 G4double G4RadioactiveDecay::GetTaoTime(const G4double t, const G4double tao)
00372 {
00373   long double taotime =0.L;
00374   G4int nbin;
00375   if ( t > SBin[NSourceBin]) {
00376     nbin  = NSourceBin;}
00377   else {
00378     nbin = 0;
00379     while (t > SBin[nbin]) nbin++;
00380     nbin--;}
00381   long double lt = t ;
00382   long double ltao = tao;
00383 
00384   if (nbin > 0) { 
00385     for (G4int i = 0; i < nbin; i++) 
00386       {
00387         taotime += (long double)SProfile[i] * (std::exp(-(lt-(long double)SBin[i+1])/ltao)-std::exp(-(lt-(long double)SBin[i])/ltao));
00388       }
00389   }
00390   taotime +=  (long double)SProfile[nbin] * (1.L-std::exp(-(lt-(long double)SBin[nbin])/ltao));
00391   if (taotime < 0.)  {
00392     G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
00393     G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
00394     G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
00395     taotime = 0.;
00396   }
00397 #ifdef G4VERBOSE
00398   if (GetVerboseLevel()>1)
00399     {G4cout <<" Tao time: " <<taotime <<G4endl;}
00400 #endif
00401   return (G4double)taotime ;
00402 }
00403 
00404 /*
00405 // Other implementation tests to avoid use of  long double
00406 G4double G4RadioactiveDecay::GetTaoTime(const G4double t, const G4double tao)
00407 {
00408   long double taotime =0.L;
00409   G4int nbin;
00410   if ( t > SBin[NSourceBin]) {
00411     nbin  = NSourceBin;}
00412   else {
00413     nbin = 0;
00414     while (t > SBin[nbin]) nbin++;
00415     nbin--;}
00416   long double lt = t ;
00417   long double ltao = tao;
00418   long double factor,factor1,dt1,dt;
00419   if (nbin > 0) {
00420             for (G4int i = 0; i < nbin; i++)
00421               { long double s1=SBin[i];
00422                 long double s2=SBin[i+1];
00423                 dt1=(s2-s1)/ltao;
00424                 if (dt1 <50.) {
00425                         factor1=std::exp(dt1)-1.;
00426                         if (factor1<dt1) factor1 =dt1;
00427                         dt=(lt-s1)/ltao;
00428                         factor=std::exp(-dt);
00429                 }
00430                 else {
00431                         factor1=1.-std::exp(-dt1);
00432                         dt=(lt-s2)/ltao;
00433                         factor=std::exp(-dt);
00434                 }
00435                 G4cout<<(long double) SProfile[i] *factor*factor1<<'\t'<<std::endl;
00436                 long double test  = (long double)SProfile[i] * (std::exp(-(lt-(long double)SBin[i+1])/ltao)-std::exp(-(lt-(long double)SBin[i])/ltao));
00437                 G4cout<<test<<std::endl;
00438                 taotime += (long double) SProfile[i] *factor*factor1;
00439            }
00440    }
00441   long double s=SBin[nbin];
00442   dt1=(lt-s)/ltao;
00443   factor=1.-std::exp(-dt1);
00444   taotime += (long double)  SProfile[nbin] *factor;
00445  if (taotime < 0.)  {
00446     G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
00447     G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
00448     G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
00449     taotime = 0.;
00450   }
00451 #ifdef G4VERBOSE
00452   if (GetVerboseLevel()>1)
00453     {G4cout <<" Tao time: " <<taotime <<G4endl;}
00454 #endif
00455   return (G4double)taotime ;
00456 }
00457 
00458 
00459 
00460 G4double G4RadioactiveDecay::GetTaoTime(const G4double t, const G4double tao)
00461 {
00462   G4double taotime =0.;
00463   G4int nbin;
00464   if ( t > SBin[NSourceBin]) {
00465     nbin  = NSourceBin;}
00466   else {
00467     nbin = 0;
00468     while (t > SBin[nbin]) nbin++;
00469     nbin--;}
00470   G4double lt = t ;
00471   G4double ltao = tao;
00472   G4double factor,factor1,dt1,dt;
00473 
00474   if (nbin > 0) {
00475     for (G4int i = 0; i < nbin; i++)
00476       { dt1=(SBin[i+1]-SBin[i])/ltao;
00477         if (dt1 <50.) {
00478                 factor1=std::exp(dt1)-1.;
00479                 if (factor1<dt1) factor1 =dt1;
00480                 dt=(lt-SBin[i])/ltao;
00481                 factor=std::exp(-(lt-SBin[i])/ltao);
00482                 G4cout<<factor<<'\t'<<factor1<<std::endl;
00483         }
00484         else {
00485                 factor1=1.-std::exp(-dt1);
00486                 factor=std::exp(-(lt-SBin[i+1])/ltao);
00487         }
00488         G4cout<<factor<<'\t'<<factor1<<std::endl;
00489         taotime += SProfile[i] *factor*factor1;
00490         G4cout<<taotime<<std::endl;
00491       }
00492   }
00493   dt1=(lt-SBin[nbin])/ltao;
00494   factor=1.-std::exp(-dt1);
00495   if (factor<(dt1-0.5*dt1*dt1)) factor =dt1-0.5*dt1*dt1;
00496 
00497 
00498   taotime +=  SProfile[nbin] *factor;
00499   G4cout<<factor<<'\t'<<taotime<<std::endl;
00500   if (taotime < 0.)  {
00501     G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
00502     G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
00503     G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
00504     taotime = 0.;
00505   }
00506 
00507 #ifdef G4VERBOSE
00508   if (GetVerboseLevel()>1)
00509     {G4cout <<" Tao time: " <<taotime <<G4endl;}
00510 #endif
00511   return (G4double)taotime ;
00512 }
00513 */
00515 //                                                                            //
00516 //  GetDecayTime                                                              //
00517 //    Randomly select a decay time for the decay process, following the       //
00518 //    supplied decay time bias scheme.                                        //
00519 //                                                                            //
00521 
00522 G4double G4RadioactiveDecay::GetDecayTime()
00523 {
00524   G4double decaytime = 0.;
00525   G4double rand = G4UniformRand();
00526   G4int i = 0;
00527   while ( DProfile[i] < rand) i++;
00528   rand = G4UniformRand();
00529   decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
00530 #ifdef G4VERBOSE
00531   if (GetVerboseLevel()>1)
00532     {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;}
00533 #endif
00534   return  decaytime;        
00535 }
00536 
00537 
00538 G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
00539 {
00540   G4int i = 0;
00541   while ( aDecayTime > DBin[i] ) i++;
00542   return  i;
00543 }
00544 
00546 //                                                                            //
00547 //  GetMeanLifeTime (required by the base class)                              //
00548 //                                                                            //
00550 
00551 G4double G4RadioactiveDecay::GetMeanLifeTime(const G4Track& theTrack,
00552                                              G4ForceCondition* )
00553 {
00554   // For varience reduction implementation the time is set to 0 so as to 
00555   // force the particle to decay immediately.
00556   // in analogueMC mode it return the particles meanlife.
00557   // 
00558   G4double meanlife = 0.;
00559   if (AnalogueMC) {
00560     const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
00561     G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
00562     G4double theLife = theParticleDef->GetPDGLifeTime();
00563 
00564 #ifdef G4VERBOSE
00565     if (GetVerboseLevel()>2)
00566       {
00567         G4cout <<"G4RadioactiveDecay::GetMeanLifeTime() " <<G4endl;
00568         G4cout <<"KineticEnergy:" <<theParticle->GetKineticEnergy()/GeV <<"[GeV]";
00569         G4cout <<"Mass:" <<theParticle->GetMass()/GeV <<"[GeV]"; 
00570         G4cout <<"Life time: " <<theLife/ns <<"[ns]" << G4endl;
00571       }
00572 #endif
00573     if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
00574     else if (theLife < 0.0) {meanlife = DBL_MAX;}
00575     else {meanlife = theLife;}
00576     // set the meanlife to zero for excited istopes which is not in the RDM database
00577     if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
00578   }
00579 #ifdef G4VERBOSE
00580   if (GetVerboseLevel()>1)
00581     {G4cout <<"mean life time: " <<meanlife/s <<"[s]" <<G4endl;}
00582 #endif
00583 
00584   return  meanlife;
00585 }
00586 
00588 //                                                                    //
00589 //  GetMeanFreePath for decay in flight                               //
00590 //                                                                    //
00592 
00593 G4double G4RadioactiveDecay::GetMeanFreePath (const G4Track& aTrack,
00594                                               G4double, G4ForceCondition*)
00595 {
00596   // get particle
00597   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00598 
00599   // returns the mean free path in GEANT4 internal units
00600   G4double pathlength;
00601   G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00602   G4double aCtau = c_light * aParticleDef->GetPDGLifeTime();
00603   G4double aMass = aParticle->GetMass();
00604 
00605 #ifdef G4VERBOSE
00606   if (GetVerboseLevel()>2) {
00607     G4cout << "G4RadioactiveDecay::GetMeanFreePath() "<< G4endl;
00608     G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
00609     G4cout << "Mass:" << aMass/GeV <<"[GeV]";
00610     G4cout << "c*Tau:" << aCtau/m <<"[m]" <<G4endl;
00611   }
00612 #endif
00613 
00614   // check if the particle is stable?
00615   if (aParticleDef->GetPDGStable()) {
00616     pathlength = DBL_MAX;
00617 
00618   } else if (aCtau < 0.0) {
00619     pathlength =  DBL_MAX;
00620 
00621     //check if the particle has very short life time ?
00622   } else if (aCtau < DBL_MIN) {
00623     pathlength =  DBL_MIN;
00624 
00625     //check if zero mass
00626   } else if (aMass <  DBL_MIN)  {
00627     pathlength =  DBL_MAX;
00628 #ifdef G4VERBOSE
00629     if (GetVerboseLevel()>1) {
00630       G4cerr << " Zero Mass particle " << G4endl;
00631     }
00632 #endif
00633   } else {
00634     //calculate the mean free path
00635     // by using normalized kinetic energy (= Ekin/mass)
00636     G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
00637     if ( rKineticEnergy > HighestValue) {
00638       // beta >> 1
00639       pathlength = ( rKineticEnergy + 1.0)* aCtau;
00640     } else if ( rKineticEnergy < DBL_MIN ) {
00641       // too slow particle
00642 #ifdef G4VERBOSE
00643       if (GetVerboseLevel()>2) {
00644         G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
00645         G4cout << aParticleDef->GetParticleName() << G4endl;
00646         G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV 
00647                <<"[GeV]";
00648       }
00649 #endif
00650       pathlength = DBL_MIN;
00651     } else {
00652       // beta << 1
00653       pathlength = aCtau*(aParticle->GetTotalMomentum())/aMass;
00654     }
00655   }
00656 #ifdef G4VERBOSE
00657   if (GetVerboseLevel()>1) {
00658     G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
00659   }
00660 #endif
00661   return  pathlength;
00662 }
00663 
00665 //                                                                    //
00666 //  BuildPhysicsTable - initialisation of atomic de-excitation        //
00667 //                                                                    //
00669 
00670 void G4RadioactiveDecay::BuildPhysicsTable(const G4ParticleDefinition&)
00671 {
00672   if(!isInitialised) {
00673     isInitialised = true;
00674     G4VAtomDeexcitation* p = G4LossTableManager::Instance()->AtomDeexcitation();
00675     if(p) { p->InitialiseAtomicDeexcitation(); }
00676   }
00677 }
00678 
00680 //                                                                            //
00681 //  LoadDecayTable loads the decay scheme from the RadioactiveDecay database  // 
00682 //  for the parent nucleus.                                                   //
00683 //                                                                            //
00685 
00686 G4DecayTable*
00687 G4RadioactiveDecay::LoadDecayTable(G4ParticleDefinition& theParentNucleus)
00688 {
00689   // Create and initialise variables used in the method.
00690   G4DecayTable* theDecayTable = new G4DecayTable();
00691 
00692   // Generate input data file name using Z and A of the parent nucleus
00693   // file containing radioactive decay data.
00694   G4int A    = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
00695   G4int Z    = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
00696   G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
00697 
00698   //Check if data have been provided by the user
00699   G4String file= theUserRadioactiveDataFiles[1000*A+Z];
00700 
00701   if (file =="") {
00702     if (!getenv("G4RADIOACTIVEDATA") ) {
00703       G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files."
00704              << G4endl;
00705       throw G4HadronicException(__FILE__, __LINE__,
00706                               "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
00707     }
00708     G4String dirName = getenv("G4RADIOACTIVEDATA");
00709 
00710     std::ostringstream os;
00711     os <<dirName <<"/z" <<Z <<".a" <<A <<'\0';
00712     file = os.str();
00713   }
00714 
00715   LoadedNuclei.push_back(theParentNucleus.GetParticleName());
00716   std::sort( LoadedNuclei.begin(), LoadedNuclei.end() );
00717   // sort needed to allow binary_search
00718 
00719   std::ifstream DecaySchemeFile(file);
00720 
00721   G4bool found(false);
00722   if (DecaySchemeFile) { 
00723     // Initialise variables used for reading in radioactive decay data.
00724     G4int nMode = 7;
00725     G4bool modeFirstRecord[7];
00726     G4double modeTotalBR[7] = {0.0};
00727     G4double modeSumBR[7];
00728     for (G4int i = 0; i < nMode; i++) {
00729       modeFirstRecord[i] = true;
00730       modeSumBR[i] = 0.0;
00731     }
00732 
00733     G4bool complete(false);
00734     char inputChars[100]={' '};
00735     G4String inputLine;
00736     G4String recordType("");
00737     G4RadioactiveDecayMode theDecayMode;
00738     G4double a(0.0);
00739     G4double b(0.0);
00740     G4double c(0.0);
00741     G4BetaDecayType betaType(allowed);
00742     G4double e0;
00743 
00744     // Loop through each data file record until you identify the decay
00745     // data relating to the nuclide of concern.
00746 
00747     while (!complete && !DecaySchemeFile.getline(inputChars, 100).eof()) {
00748       inputLine = inputChars;
00749       inputLine = inputLine.strip(1);
00750       if (inputChars[0] != '#' && inputLine.length() != 0) {
00751         std::istringstream tmpStream(inputLine);
00752 
00753         if (inputChars[0] == 'P') {
00754           // Nucleus is a parent type.  Check excitation level to see if it
00755           // matches that of theParentNucleus
00756           tmpStream >> recordType >> a >> b;
00757           if (found) {complete = true;}
00758           else {found = (std::abs(a*keV - E) < levelTolerance);}
00759 
00760         } else if (found) {
00761           // The right part of the radioactive decay data file has been found.  Search
00762           // through it to determine the mode of decay of the subsequent records.
00763           if (inputChars[0] == 'W') {
00764 #ifdef G4VERBOSE
00765             if (GetVerboseLevel() > 0) {
00766               // a comment line identified and print out the message
00767               //
00768               G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
00769               G4cout << "   In data file " << file << G4endl;
00770               G4cout << "   " << inputLine << G4endl;
00771             }
00772 #endif
00773           } else {
00774             tmpStream >> theDecayMode >> a >> b >> c >> betaType;
00775 
00776             // Allowed transitions are the default. Forbidden transitions are
00777             // indicated in the last column.
00778             if (inputLine.length() < 80) betaType = allowed;
00779             a /= 1000.;
00780             c /= 1000.;
00781 
00782             switch (theDecayMode) {
00783 
00784               case IT:   // Isomeric transition
00785               {
00786                 G4ITDecayChannel* anITChannel =
00787                   new G4ITDecayChannel(GetVerboseLevel(),
00788                                        (const G4Ions*)& theParentNucleus, b);
00789                 anITChannel->SetICM(applyICM);
00790                 anITChannel->SetARM(applyARM);
00791                 anITChannel->SetHLThreshold(halflifethreshold);
00792                 theDecayTable->Insert(anITChannel);
00793               }
00794               break;
00795 
00796               case BetaMinus:
00797               {
00798                 if (modeFirstRecord[1]) {
00799                   modeFirstRecord[1] = false;
00800                   modeTotalBR[1] = b;
00801                 } else {
00802                   if (c > 0.) {
00803                     e0 = c*MeV/0.511;
00804                     G4BetaDecayCorrections corrections(Z+1, A);
00805 
00806                     // array to store function shape
00807                     G4int npti = 100;                           
00808                     G4double* pdf = new G4double[npti];
00809 
00810                     G4double e;   // Total electron energy in units of electron mass
00811                     G4double p;   // Electron momentum in units of electron mass 
00812                     G4double f;   // Spectral shape function value
00813                     for (G4int ptn = 0; ptn < npti; ptn++) {
00814                       // Calculate simple phase space spectrum
00815                       e = 1. + e0*(ptn+0.5)/100.;
00816                       p = std::sqrt(e*e - 1.);
00817                       f = p*e*(e0-e+1)*(e0-e+1);
00818 
00819                       // Apply Fermi factor to get allowed shape
00820                       f *= corrections.FermiFunction(e);
00821 
00822                       // Apply shape factor for forbidden transitions
00823                       f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
00824                       pdf[ptn] = f;
00825                     }
00826 
00827                     RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
00828                     G4BetaMinusDecayChannel *aBetaMinusChannel = new
00829                     G4BetaMinusDecayChannel(GetVerboseLevel(), &theParentNucleus,
00830                                             b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
00831                     aBetaMinusChannel->SetICM(applyICM);
00832                     aBetaMinusChannel->SetARM(applyARM);
00833                     aBetaMinusChannel->SetHLThreshold(halflifethreshold);
00834                     theDecayTable->Insert(aBetaMinusChannel);
00835                     modeSumBR[1] += b;
00836                     delete[] pdf;
00837                   } // c > 0
00838                 } // if not first record
00839               }
00840               break;
00841 
00842               case BetaPlus:
00843               {
00844                 if (modeFirstRecord[2]) {
00845                   modeFirstRecord[2] = false;
00846                   modeTotalBR[2] = b;
00847                 } else {
00848                   e0 = c*MeV/0.511 - 2.;
00849                   // Need to test e0 for nuclei which have Q < 2Me in their
00850                   // data files (e.g. z67.a162)
00851                   if (e0 > 0.) {
00852                     G4BetaDecayCorrections corrections(1-Z, A);
00853 
00854                     // array to store function shape
00855                     G4int npti = 100;                           
00856                     G4double* pdf = new G4double[npti];
00857 
00858                     G4double e;   // Total positron energy in units of electron mass
00859                     G4double p;   // Positron momentum in units of electron mass
00860                     G4double f;   // Spectral shape function value
00861                     for (G4int ptn = 0; ptn < npti; ptn++) {
00862                       // Calculate simple phase space spectrum
00863                       e = 1. + e0*(ptn+0.5)/100.;
00864                       p = std::sqrt(e*e - 1.);
00865                       f = p*e*(e0-e+1)*(e0-e+1);
00866 
00867                       // Apply Fermi factor to get allowed shape
00868                       f *= corrections.FermiFunction(e);
00869 
00870                       // Apply shape factor for forbidden transitions
00871                       f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
00872                       pdf[ptn] = f;
00873                     }
00874                     RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
00875                     G4BetaPlusDecayChannel *aBetaPlusChannel = new 
00876                     G4BetaPlusDecayChannel(GetVerboseLevel(), &theParentNucleus,
00877                                            b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
00878                     aBetaPlusChannel->SetICM(applyICM);
00879                     aBetaPlusChannel->SetARM(applyARM);
00880                     aBetaPlusChannel->SetHLThreshold(halflifethreshold);
00881                     theDecayTable->Insert(aBetaPlusChannel);
00882                     modeSumBR[2] += b;
00883                     delete[] pdf;
00884                   } // if e0 > 0
00885                 } // if not first record
00886               }
00887               break;
00888 
00889               case KshellEC:  // K-shell electron capture
00890 
00891                 if (modeFirstRecord[3]) {
00892                   modeFirstRecord[3] = false;
00893                   modeTotalBR[3] = b;
00894                 } else {
00895                   G4KshellECDecayChannel* aKECChannel =
00896                     new G4KshellECDecayChannel(GetVerboseLevel(),
00897                                                &theParentNucleus,
00898                                                b, c*MeV, a*MeV);
00899                   aKECChannel->SetICM(applyICM);
00900                   aKECChannel->SetARM(applyARM);
00901                   aKECChannel->SetHLThreshold(halflifethreshold);
00902                   theDecayTable->Insert(aKECChannel);
00903                   modeSumBR[3] += b;
00904                 }
00905                 break;
00906 
00907               case LshellEC:  // L-shell electron capture
00908 
00909                 if (modeFirstRecord[4]) {
00910                   modeFirstRecord[4] = false;
00911                   modeTotalBR[4] = b;
00912                 } else {
00913                   G4LshellECDecayChannel *aLECChannel = new
00914                     G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
00915                                             b, c*MeV, a*MeV);
00916                   aLECChannel->SetICM(applyICM);
00917                   aLECChannel->SetARM(applyARM);
00918                   aLECChannel->SetHLThreshold(halflifethreshold);
00919                   theDecayTable->Insert(aLECChannel);
00920                   modeSumBR[4] += b;
00921                 }
00922                 break;
00923 
00924               case MshellEC:  // M-shell electron capture
00925                               // In this implementation it is added to L-shell case
00926                 if (modeFirstRecord[5]) {
00927                   modeFirstRecord[5] = false;
00928                   modeTotalBR[5] = b;
00929                 } else {
00930                   G4MshellECDecayChannel* aMECChannel =
00931                     new G4MshellECDecayChannel(GetVerboseLevel(),
00932                                                &theParentNucleus,
00933                                                b, c*MeV, a*MeV);
00934                   aMECChannel->SetICM(applyICM);
00935                   aMECChannel->SetARM(applyARM);
00936                   aMECChannel->SetHLThreshold(halflifethreshold);
00937                   theDecayTable->Insert(aMECChannel);
00938                   modeSumBR[5] += b;
00939                 }
00940                 break;
00941 
00942               case Alpha:
00943                   //G4cout<<"Alpha channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
00944 
00945                   if (modeFirstRecord[6]) {
00946                   modeFirstRecord[6] = false;
00947                   modeTotalBR[6] = b;
00948                 } else {
00949                   G4AlphaDecayChannel* anAlphaChannel =
00950                      new G4AlphaDecayChannel(GetVerboseLevel(),
00951                                              &theParentNucleus,
00952                                              b, c*MeV, a*MeV);
00953                   anAlphaChannel->SetICM(applyICM);
00954                   anAlphaChannel->SetARM(applyARM);
00955                   anAlphaChannel->SetHLThreshold(halflifethreshold);
00956                   theDecayTable->Insert(anAlphaChannel);
00957                   modeSumBR[6] += b;
00958                 }
00959                 break;
00960               case SpFission:
00961                   //Still needed to be implemented
00962                   //G4cout<<"Sp fission channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
00963                   break;
00964               case ERROR:
00965 
00966               default:
00967                 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
00968                             FatalException, "Selected decay mode does not exist");
00969             }  // switch
00970           }  // if char == W
00971         }  // if char == P 
00972       }  // if char != #
00973     }  // While
00974 
00975     // Go through the decay table and make sure that the branching ratios are
00976     // correctly normalised.
00977 
00978     G4VDecayChannel* theChannel = 0;
00979     G4NuclearDecayChannel* theNuclearDecayChannel = 0;
00980     G4String mode = "";
00981 
00982     G4double theBR = 0.0;
00983     for (G4int i = 0; i < theDecayTable->entries(); i++) {
00984       theChannel = theDecayTable->GetDecayChannel(i);
00985       theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
00986       theDecayMode = theNuclearDecayChannel->GetDecayMode();
00987 
00988       if (theDecayMode != IT) {
00989         theBR = theChannel->GetBR();
00990         theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
00991       }
00992     } 
00993   }   // if (DecaySchemeFile)   
00994   DecaySchemeFile.close();
00995 
00996                 
00997   if (!found && E > 0.) {
00998     // cases where IT cascade for exited isotopes without entry in RDM database
00999     // Decay mode is isomeric transition.
01000     //
01001     G4ITDecayChannel *anITChannel = new G4ITDecayChannel
01002       (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
01003     anITChannel->SetICM(applyICM);
01004     anITChannel->SetARM(applyARM);
01005     anITChannel->SetHLThreshold(halflifethreshold);
01006     theDecayTable->Insert(anITChannel);
01007   } 
01008   if (!theDecayTable) {
01009     //
01010     // There is no radioactive decay data for this nucleus.  Return a null
01011     // decay table.
01012     //
01013     G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
01014     theDecayTable = 0;
01015     return theDecayTable;
01016   }     
01017   if (theDecayTable && GetVerboseLevel()>1)
01018     {
01019       G4cout <<"G4RadioactiveDecay::LoadDecayTable()" << G4endl;
01020       G4cout << "  No. of  entries: "<< theDecayTable->entries() <<G4endl;
01021       theDecayTable ->DumpInfo();
01022     }
01023 
01024   return theDecayTable;
01025 }
01027 //
01028 void G4RadioactiveDecay::AddUserDecayDataFile(G4int Z, G4int A,G4String filename)
01029 { if (Z<1 || A<2) {
01030         G4cout<<"Z and A not valid!"<<G4endl;
01031   }
01032 
01033   std::ifstream DecaySchemeFile(filename);
01034   if (DecaySchemeFile){
01035         G4int ID_ion=A*1000+Z;
01036         theUserRadioactiveDataFiles[ID_ion]=filename;
01037         theIsotopeTable->AddUserDecayDataFile(Z,A,filename);
01038   }
01039   else {
01040         G4cout<<"The file "<<filename<<" does not exist!"<<G4endl;
01041   }
01042 }
01044 //
01045 
01046 
01047 void
01048 G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE, 
01049                                  G4int theG, std::vector<G4double> theRates, 
01050                                  std::vector<G4double> theTaos)
01051 { 
01052   //fill the decay rate vector 
01053   theDecayRate.SetZ(theZ);
01054   theDecayRate.SetA(theA);
01055   theDecayRate.SetE(theE);
01056   theDecayRate.SetGeneration(theG);
01057   theDecayRate.SetDecayRateC(theRates);
01058   theDecayRate.SetTaos(theTaos);
01059 }
01060 
01062 // 
01063 void
01064 G4RadioactiveDecay::AddDecayRateTable(const G4ParticleDefinition& theParentNucleus)
01065 {
01066   // 1) To calculate all the coefficiecies required to derive the
01067   //    radioactivities for all progeny of theParentNucleus
01068   //
01069   // 2) Add the coefficiencies to the decay rate table vector 
01070   //
01071 
01072   //
01073   // Create and initialise variables used in the method.
01074   //
01075   theDecayRateVector.clear();
01076 
01077   G4int nGeneration = 0;
01078   std::vector<G4double> rates;
01079   std::vector<G4double> taos;
01080 
01081   // start rate is -1.
01082   // Eq.4.26 of the Technical Note
01083   rates.push_back(-1.);
01084   //
01085   //
01086   G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
01087   G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
01088   G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
01089   G4double tao = theParentNucleus.GetPDGLifeTime();
01090   if (tao < 0.) tao = 1e-100;
01091   taos.push_back(tao);
01092   G4int nEntry = 0;
01093 
01094   //fill the decay rate with the intial isotope data
01095   SetDecayRate(Z,A,E,nGeneration,rates,taos);
01096 
01097   // store the decay rate in decay rate vector
01098   theDecayRateVector.push_back(theDecayRate);
01099   nEntry++;
01100 
01101   // now start treating the sencondary generations..
01102 
01103   G4bool stable = false;
01104   G4int i;
01105   G4int j;
01106   G4VDecayChannel* theChannel = 0;
01107   G4NuclearDecayChannel* theNuclearDecayChannel = 0;
01108   G4ITDecayChannel* theITChannel = 0;
01109   G4BetaMinusDecayChannel *theBetaMinusChannel = 0;
01110   G4BetaPlusDecayChannel *theBetaPlusChannel = 0;
01111   G4AlphaDecayChannel *theAlphaChannel = 0;
01112   G4RadioactiveDecayMode theDecayMode;
01113   G4double theBR = 0.0;
01114   G4int AP = 0;
01115   G4int ZP = 0;
01116   G4int AD = 0;
01117   G4int ZD = 0;
01118   G4double EP = 0.;
01119   std::vector<G4double> TP;
01120   std::vector<G4double> RP;
01121   G4ParticleDefinition *theDaughterNucleus;
01122   G4double daughterExcitation;
01123   G4ParticleDefinition *aParentNucleus;
01124   G4IonTable* theIonTable;
01125   G4DecayTable *aTempDecayTable;
01126   G4double theRate;
01127   G4double TaoPlus;
01128   G4int nS = 0;
01129   G4int nT = nEntry;
01130   G4double brs[7];
01131   //
01132   theIonTable =
01133     (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
01134 
01135   while (!stable) {
01136     nGeneration++;
01137     for (j = nS; j< nT; j++) {
01138       ZP = theDecayRateVector[j].GetZ();
01139       AP = theDecayRateVector[j].GetA();
01140       EP = theDecayRateVector[j].GetE();
01141       RP = theDecayRateVector[j].GetDecayRateC();
01142       TP = theDecayRateVector[j].GetTaos();      
01143       if (GetVerboseLevel()>0){
01144         G4cout <<"G4RadioactiveDecay::AddDecayRateTable : "
01145                << " daughters of ("<< ZP <<", "<<AP<<", "
01146                << EP <<") "
01147                << " are being calculated. "       
01148                <<" generation = "
01149                << nGeneration << G4endl;
01150       }
01151       aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
01152       if (!IsLoaded(*aParentNucleus)){
01153         aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
01154       }
01155         
01156       G4DecayTable* theDecayTable = new G4DecayTable();
01157       aTempDecayTable = aParentNucleus->GetDecayTable();
01158       for (i=0; i< 7; i++) brs[i] = 0.0;
01159 
01160       //
01161       // Go through the decay table and to combine the same decay channels
01162       //
01163       for (i=0; i<aTempDecayTable->entries(); i++) {
01164         theChannel             = aTempDecayTable->GetDecayChannel(i);
01165         theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
01166         theDecayMode           = theNuclearDecayChannel->GetDecayMode();
01167         daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation ();
01168         theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus () ;
01169         AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
01170         ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();  
01171         G4NuclearLevelManager* levelManager = 
01172            G4NuclearLevelStore::GetInstance()->GetManager(ZD, AD);
01173         if (levelManager->NumberOfLevels() ) {
01174           const G4NuclearLevel* level =
01175               levelManager->NearestLevel (daughterExcitation);
01176 
01177           if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
01178             // Level half-life is in ns and the threshold is set to 1 micros by default, user can set it via the UI command
01179             if (level->HalfLife()*ns >= halflifethreshold ){    
01180               // save the metastable nucleus 
01181               theDecayTable->Insert(theChannel);
01182             } 
01183             else{
01184               brs[theDecayMode] += theChannel->GetBR();
01185             }
01186           }
01187           else {
01188             brs[theDecayMode] += theChannel->GetBR();
01189           }
01190         }
01191         else{
01192           brs[theDecayMode] += theChannel->GetBR();
01193         }
01194       }     
01195       brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
01196       brs[3] = brs[4] =brs[5] =  0.0;
01197       for (i= 0; i<7; i++){
01198         if (brs[i] > 0.) {
01199           switch ( i ) {
01200           case 0:
01201             //
01202             //
01203             // Decay mode is isomeric transition.
01204             //
01205 
01206             theITChannel =  new G4ITDecayChannel
01207               (0, (const G4Ions*) aParentNucleus, brs[0]);
01208 
01209             theDecayTable->Insert(theITChannel);
01210             break;
01211 
01212           case 1:
01213             //
01214             //
01215             // Decay mode is beta-.
01216             //
01217             theBetaMinusChannel = new G4BetaMinusDecayChannel (0, aParentNucleus,
01218                                                                brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
01219             theDecayTable->Insert(theBetaMinusChannel);
01220 
01221             break;
01222 
01223           case 2:
01224             //
01225             //
01226             // Decay mode is beta+ + EC.
01227             //
01228             theBetaPlusChannel = new G4BetaPlusDecayChannel (GetVerboseLevel(), aParentNucleus,
01229                                                              brs[2], 0.*MeV, 0.*MeV, 1, false, 0);
01230             theDecayTable->Insert(theBetaPlusChannel);
01231             break;                    
01232 
01233           case 6:
01234             //
01235             //
01236             // Decay mode is alpha.
01237             //
01238             theAlphaChannel = new G4AlphaDecayChannel(GetVerboseLevel(),
01239                                                       aParentNucleus,
01240                                                       brs[6], 0.*MeV, 0.*MeV);
01241             theDecayTable->Insert(theAlphaChannel);
01242             break;
01243 
01244           default:
01245             break;
01246           }
01247         }
01248       }
01249       // 
01250       // loop over all branches in theDecayTable
01251       //
01252       for (i = 0; i < theDecayTable->entries(); i++){
01253         theChannel = theDecayTable->GetDecayChannel(i);
01254         theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
01255         theBR = theChannel->GetBR();
01256         theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
01257         //  first check if the decay of the original nucleus is an IT channel, if true create a new groud-level nucleus
01258         if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
01259           A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
01260           Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
01261           theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
01262         }
01263         if (IsApplicable(*theDaughterNucleus) &&
01264             theBR && 
01265             aParentNucleus != theDaughterNucleus) { 
01266           // need to make sure daugher has decaytable
01267           if (!IsLoaded(*theDaughterNucleus)){
01268             theDaughterNucleus->SetDecayTable(LoadDecayTable(*theDaughterNucleus));
01269           }
01270           if (theDaughterNucleus->GetDecayTable()->entries() ) {
01271             //
01272             A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
01273             Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
01274             E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
01275 
01276             TaoPlus = theDaughterNucleus->GetPDGLifeTime();
01277             //          cout << TaoPlus <<G4endl;
01278 
01279             if (TaoPlus <= 0.)  TaoPlus = 1e-100;
01280 
01281 
01282 
01283             // first set the taos, one simply need to add to the parent ones
01284             taos.clear();
01285             taos = TP;
01286             size_t k;
01287             //check that TaoPlus differs from other taos from at least 1.e5 relative difference
01288             //for (k = 0; k < TP.size(); k++){
01289             //  if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
01290             //}
01291             taos.push_back(TaoPlus);
01292             // now calculate the coefficiencies
01293             //
01294             // they are in two parts, first the less than n ones
01295             // Eq 4.24 of the TN
01296             rates.clear();
01297             long double ta1,ta2;
01298             ta2 = (long double)TaoPlus;
01299             for (k = 0; k < RP.size(); k++){
01300               ta1 = (long double)TP[k];
01301               if (ta1 == ta2) {
01302                 theRate = 1.e100;
01303               }else{
01304                 theRate = ta1/(ta1-ta2);}
01305               theRate = theRate * theBR * RP[k];
01306               rates.push_back(theRate);
01307             }
01308             //
01309             // the sencond part: the n:n coefficiency
01310             // Eq 4.25 of the TN.  Note Yn+1 is zero apart from Y1 which is -1 as treated at line 1013
01311             // 
01312             theRate = 0.;
01313             long double aRate, aRate1;
01314             aRate1 = 0.L;
01315             for (k = 0; k < RP.size(); k++){
01316               ta1 = (long double)TP[k];
01317               if (ta1 == ta2 ) {
01318                 aRate = 1.e100;
01319               }else {
01320                 aRate = ta2/(ta1-ta2);}
01321               aRate = aRate * (long double)(theBR * RP[k]);
01322               aRate1 += aRate;
01323             }
01324             theRate = -aRate1;
01325             rates.push_back(theRate);         
01326             SetDecayRate (Z,A,E,nGeneration,rates,taos);
01327             theDecayRateVector.push_back(theDecayRate);
01328             nEntry++;
01329           }
01330         } // end of testing daughter nucleus
01331       } // end of i loop( the branches) 
01332       //      delete theDecayTable;
01333 
01334     } //end of for j loop
01335     nS = nT;
01336     nT = nEntry;
01337     if (nS == nT) stable = true;
01338   }
01339 
01340   //end of while loop
01341   // the calculation completed here
01342 
01343 
01344   // fill the first part of the decay rate table
01345   // which is the name of the original particle (isotope) 
01346   //
01347   theDecayRateTable.SetIonName(theParentNucleus.GetParticleName()); 
01348   //
01349   //
01350   // now fill the decay table with the newly completed decay rate vector
01351   //
01352 
01353   theDecayRateTable.SetItsRates(theDecayRateVector);
01354 
01355   //
01356   // finally add the decayratetable to the tablevector
01357   //
01358   theDecayRateTableVector.push_back(theDecayRateTable);
01359 }
01360 
01362 //                                                                            //
01363 //  SetSourceTimeProfile                                                      //
01364 //     read in the source time profile function (histogram)                   //
01365 //                                                                            //
01367 
01368 void G4RadioactiveDecay::SetSourceTimeProfile(G4String filename)
01369 {
01370   std::ifstream infile ( filename, std::ios::in );
01371   if (!infile) {
01372     G4ExceptionDescription ed;
01373     ed << " Could not open file " << filename << G4endl; 
01374     G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
01375                 FatalException, ed);
01376   }
01377 
01378   G4double bin, flux;
01379   NSourceBin = -1;
01380   while (infile >> bin >> flux ) {
01381     NSourceBin++;
01382     if (NSourceBin > 99) {
01383       G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
01384                   FatalException, "Input source time file too big (>100 rows)");
01385 
01386     } else {
01387       SBin[NSourceBin] = bin * s;
01388       SProfile[NSourceBin] = flux;
01389     }
01390   }
01391   SetAnalogueMonteCarlo(0);
01392   infile.close();
01393 
01394 #ifdef G4VERBOSE
01395   if (GetVerboseLevel()>1)
01396     {G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;}
01397 #endif
01398 }
01399 
01401 //                                                                            //
01402 //  SetDecayBiasProfile                                                       //
01403 //    read in the decay bias scheme function (histogram)                      //
01404 //                                                                            //
01406 
01407 void G4RadioactiveDecay::SetDecayBias(G4String filename)
01408 {
01409   
01410   std::ifstream infile(filename, std::ios::in);
01411   if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_003",
01412                            FatalException, "Unable to open bias data file" );
01413 
01414   G4double bin, flux;
01415   G4int dWindows = 0;
01416   G4int i ;
01417 
01418   theRadioactivityTables.clear();
01419 
01420   NDecayBin = -1;
01421   while (infile >> bin >> flux ) {
01422     NDecayBin++;
01423     if (NDecayBin > 99) {
01424       G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_004",
01425                   FatalException, "Input bias file too big (>100 rows)" );
01426     } else {
01427       DBin[NDecayBin] = bin * s;
01428       DProfile[NDecayBin] = flux;
01429       if (flux > 0.) {
01430         decayWindows[NDecayBin] = dWindows;
01431         dWindows++;
01432         G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
01433         theRadioactivityTables.push_back(rTable);
01434       }
01435     }
01436   }
01437   for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
01438   for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
01439   // converted to accumulated probabilities
01440 
01441   SetAnalogueMonteCarlo(0);
01442   infile.close();
01443 
01444 #ifdef G4VERBOSE
01445   if (GetVerboseLevel()>1)
01446     {G4cout <<" Decay Bias Profile  Nbin = " << NDecayBin <<G4endl;}
01447 #endif
01448 }
01449 
01451 //                                                                            //
01452 //  DecayIt                                                                   //
01453 //                                                                            //
01455 
01456 G4VParticleChange*
01457 G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
01458 {
01459   // Initialize the G4ParticleChange object. Get the particle details and the
01460   // decay table.
01461 
01462 
01463 
01464   fParticleChangeForRadDecay.Initialize(theTrack);
01465   const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
01466 
01467 
01468   G4ParticleDefinition *theParticleDef = theParticle->GetDefinition();
01469 
01470 
01471   // First check whether RDM applies to the current logical volume
01472   if (!isAllVolumesMode){
01473    if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
01474                           theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
01475 #ifdef G4VERBOSE
01476     if (GetVerboseLevel()>0) {
01477       G4cout <<"G4RadioactiveDecay::DecayIt : "
01478              << theTrack.GetVolume()->GetLogicalVolume()->GetName()
01479              << " is not selected for the RDM"<< G4endl;
01480       G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
01481       G4cout << " The Valid volumes are " << G4endl;
01482       for (size_t i = 0; i< ValidVolumes.size(); i++) G4cout << ValidVolumes[i] << G4endl;
01483     }
01484 #endif
01485     fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01486 
01487     // Kill the parent particle.
01488 
01489     fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01490     fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01491     ClearNumberOfInteractionLengthLeft();
01492     return &fParticleChangeForRadDecay;
01493    }
01494   }
01495 
01496   // now check is the particle is valid for RDM
01497 
01498   if (!(IsApplicable(*theParticleDef))) { 
01499     //
01500     // The particle is not a Ion or outside the nucleuslimits for decay
01501     //
01502 #ifdef G4VERBOSE
01503     if (GetVerboseLevel()>0) {
01504       G4cerr <<"G4RadioactiveDecay::DecayIt : "
01505              <<theParticleDef->GetParticleName() 
01506              << " is not a valid nucleus for the RDM"<< G4endl;
01507     }
01508 #endif
01509     fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01510 
01511     //
01512     // Kill the parent particle.
01513     //
01514     fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01515     fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01516     ClearNumberOfInteractionLengthLeft();
01517     return &fParticleChangeForRadDecay;
01518   }
01519 
01520   if (!IsLoaded(*theParticleDef))
01521     theParticleDef->SetDecayTable(LoadDecayTable(*theParticleDef));
01522 
01523   G4DecayTable* theDecayTable = theParticleDef->GetDecayTable();
01524 
01525   if (theDecayTable == 0 || theDecayTable->entries() == 0) {
01526       // There are no data in the decay table.  Set the particle change parameters
01527       // to indicate this.
01528 #ifdef G4VERBOSE
01529       if (GetVerboseLevel()>0)
01530         {
01531           G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined  for ";
01532           G4cerr <<theParticleDef->GetParticleName() <<G4endl;
01533         }
01534 #endif
01535       fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01536 
01537       // Kill the parent particle.
01538       fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01539       fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01540       ClearNumberOfInteractionLengthLeft();
01541       return &fParticleChangeForRadDecay;
01542 
01543   } else { 
01544     // Data found.  Try to decay nucleus
01545     G4double energyDeposit = 0.0;
01546     G4double finalGlobalTime = theTrack.GetGlobalTime();
01547     G4double finalLocalTime = theTrack.GetLocalTime();
01548     G4int index;
01549     G4ThreeVector currentPosition;
01550     currentPosition = theTrack.GetPosition();
01551 
01552     // Check whether use Analogue or VR implementation
01553     if (AnalogueMC) {
01554 #ifdef G4VERBOSE
01555       if (GetVerboseLevel() > 0)
01556         G4cout <<"DecayIt:  Analogue MC version " << G4endl;
01557 #endif
01558 
01559       G4DecayProducts* products = DoDecay(*theParticleDef);
01560 
01561       // Check if the product is the same as input and kill the track if
01562       // necessary to prevent infinite loop (11/05/10, F.Lei)
01563       if ( products->entries() == 1) {
01564         fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01565         fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01566         fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01567         ClearNumberOfInteractionLengthLeft();
01568         return &fParticleChangeForRadDecay;
01569       }
01570 
01571       // Get parent particle information and boost the decay products to the
01572       // laboratory frame based on this information.
01573 
01574 
01575       //The Parent Energy used for the boost should be the total energy of
01576       // the nucleus of the parent ion without the energy of the shell electrons
01577       // (correction for bug 1359 by L. Desorgher)
01578       G4double ParentEnergy =  theParticle->GetKineticEnergy()+theParticle->GetParticleDefinition()->GetPDGMass();
01579       G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
01580 
01581 
01582 
01583       if (theTrack.GetTrackStatus() == fStopButAlive) { //this condition seems to be always True, further investigation is needed (L.Desorgher)
01584 
01585         // The particle is decayed at rest.
01586         // since the time is still for rest particle in G4 we need to add the
01587         // additional time lapsed between the particle come to rest and the
01588         // actual decay.  This time is simply sampled with the mean-life of
01589         // the particle.  But we need to protect the case PDGTime < 0.
01590         // (F.Lei 11/05/10)
01591         G4double temptime = -std::log( G4UniformRand())
01592                             *theParticleDef->GetPDGLifeTime();
01593         if (temptime < 0.) temptime = 0.; 
01594         finalGlobalTime += temptime;
01595         finalLocalTime += temptime;
01596         energyDeposit += theParticle->GetKineticEnergy();
01597       }
01598       products->Boost( ParentEnergy, ParentDirection);
01599 
01600 
01601 
01602       // Add products in theParticleChangeForRadDecay.
01603       G4int numberOfSecondaries = products->entries();
01604       fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
01605 #ifdef G4VERBOSE
01606       if (GetVerboseLevel()>1) {
01607         G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
01608         G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
01609         G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
01610         G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
01611         G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
01612         G4cout << G4endl;
01613         G4cout <<"G4Decay::DecayIt  : decay products in Lab. Frame" <<G4endl;
01614         products->DumpInfo();
01615         products->IsChecked();
01616       }
01617 #endif
01618       for (index=0; index < numberOfSecondaries; index++) {
01619         G4Track* secondary = new G4Track(products->PopProducts(),
01620                                          finalGlobalTime, currentPosition);
01621         secondary->SetGoodForTrackingFlag();
01622         secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
01623         fParticleChangeForRadDecay.AddSecondary(secondary);
01624       }
01625       delete products;
01626       // end of analogue MC algorithm
01627 
01628     } else {
01629       // Variance Reduction Method
01630 #ifdef G4VERBOSE
01631       if (GetVerboseLevel()>0)
01632         G4cout << "DecayIt: Variance Reduction version " << G4endl;
01633 #endif
01634       if (!IsRateTableReady(*theParticleDef)) {
01635         // if the decayrates are not ready, calculate them and 
01636         // add to the rate table vector 
01637         AddDecayRateTable(*theParticleDef);
01638       }
01639       //retrieve the rates 
01640       GetDecayRateTable(*theParticleDef);
01641 
01642       // declare some of the variables required in the implementation
01643       G4ParticleDefinition* parentNucleus;
01644       G4IonTable* theIonTable;
01645       G4int PZ;
01646       G4int PA;
01647       G4double PE;
01648       std::vector<G4double> PT;
01649       std::vector<G4double> PR;
01650       G4double taotime;
01651       long double decayRate;
01652 
01653       size_t i;
01654       size_t j;
01655       G4int numberOfSecondaries;
01656       G4int totalNumberOfSecondaries = 0;
01657       G4double currentTime = 0.;
01658       G4int ndecaych;
01659       G4DynamicParticle* asecondaryparticle;
01660       std::vector<G4DynamicParticle*> secondaryparticles;
01661       std::vector<G4double> pw;
01662       std::vector<G4double> ptime;
01663       pw.clear();
01664       ptime.clear();
01665 
01666       //now apply the nucleus splitting
01667       for (G4int n = 0; n < NSplit; n++) {
01668         // Get the decay time following the decay probability function 
01669         // suppllied by user  
01670         G4double theDecayTime = GetDecayTime();
01671         G4int nbin = GetDecayTimeBin(theDecayTime);
01672             
01673         // calculate the first part of the weight function  
01674         G4double weight1 = 1.; 
01675         if (nbin == 1) {
01676           weight1 = 1./DProfile[nbin-1] 
01677                     *(DBin[nbin]-DBin[nbin-1])/NSplit;
01678         } else if (nbin > 1) {
01679           weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
01680                     *(DBin[nbin]-DBin[nbin-1])/NSplit;
01681         }
01682 
01683         // it should be calculated in seconds
01684         weight1 /= s ;
01685             
01686         // loop over all the possible secondaries of the nucleus
01687         // the first one is itself.
01688         for (i = 0; i < theDecayRateVector.size(); i++) {
01689           PZ = theDecayRateVector[i].GetZ();
01690           PA = theDecayRateVector[i].GetA();
01691           PE = theDecayRateVector[i].GetE();
01692           PT = theDecayRateVector[i].GetTaos();
01693           PR = theDecayRateVector[i].GetDecayRateC();
01694 
01695           // Calculate the decay rate of the isotope
01696           // decayRate is the radioactivity of isotope (PZ,PA,PE) at the 
01697           // time 'theDecayTime'
01698           // it will be used to calculate the statistical weight of the 
01699           // decay products of this isotope
01700 
01701           //  G4cout <<"PA= "<< PA << " PZ= " << PZ << " PE= "<< PE  <<G4endl;
01702           decayRate = 0.L;
01703           for (j = 0; j < PT.size(); j++) {
01704             taotime = GetTaoTime(theDecayTime,PT[j]);
01705             decayRate -= PR[j] * (long double)taotime;
01706             // Eq.4.23 of of the TN
01707             // note the negative here is required as the rate in the
01708             // equation is defined to be negative, 
01709             // i.e. decay away, but we need positive value here.
01710 
01711             // G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t"
01712             //        << decayRate << G4endl;           
01713           }
01714 
01715           // add the isotope to the radioactivity tables
01716           //  G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
01717           //  G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
01718           theRadioactivityTables[decayWindows[nbin-1]]->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
01719 
01720           // Now calculate the statistical weight
01721           // One needs to fold the source bias function with the decaytime
01722           // also need to include the track weight! (F.Lei, 28/10/10)
01723           G4double weight = weight1*decayRate*theTrack.GetWeight();
01724 
01725 
01726 
01727           // decay the isotope 
01728           theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
01729           parentNucleus = theIonTable->GetIon(PZ,PA,PE);
01730 
01731           // Create a temprary products buffer.
01732           // Its contents to be transfered to the products at the end of the loop
01733           G4DecayProducts* tempprods = 0;
01734 
01735           // Decide whether to apply branching ratio bias or not             
01736           if (BRBias) {
01737             G4DecayTable* decayTable = parentNucleus->GetDecayTable();
01738             ndecaych = G4int(decayTable->entries()*G4UniformRand());
01739             G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
01740             if (theDecayChannel == 0) {
01741               // Decay channel not found.
01742 #ifdef G4VERBOSE
01743               if (GetVerboseLevel()>0) {
01744                 G4cerr << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
01745                 G4cerr << " for this nucleus; decay as if no biasing active ";
01746                 G4cerr << G4endl;
01747                 decayTable ->DumpInfo();
01748               }
01749 #endif
01750               tempprods = DoDecay(*parentNucleus);  // DHW 6 Dec 2010 - do decay as if no biasing
01751                                                     //           to avoid deref of temppprods = 0
01752             } else {
01753               // A decay channel has been identified, so execute the DecayIt.
01754               G4double tempmass = parentNucleus->GetPDGMass();
01755               tempprods = theDecayChannel->DecayIt(tempmass);
01756               weight *= (theDecayChannel->GetBR())*(decayTable->entries());
01757             }
01758           } else {
01759             tempprods = DoDecay(*parentNucleus);
01760           }
01761 
01762           // save the secondaries for buffers
01763           numberOfSecondaries = tempprods->entries();
01764           currentTime = finalGlobalTime + theDecayTime;
01765           for (index = 0; index < numberOfSecondaries; index++) {
01766             asecondaryparticle = tempprods->PopProducts();
01767             if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5) {
01768               pw.push_back(weight);
01769               ptime.push_back(currentTime);
01770               secondaryparticles.push_back(asecondaryparticle);
01771             }
01772           }
01773           delete tempprods;
01774 
01775         } // end of i loop
01776       } // end of n loop 
01777 
01778       // now deal with the secondaries in the two stl containers
01779       // and submmit them back to the tracking manager
01780       totalNumberOfSecondaries = pw.size();
01781       fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
01782       for (index=0; index < totalNumberOfSecondaries; index++) { 
01783         G4Track* secondary = new G4Track(secondaryparticles[index],
01784                                          ptime[index], currentPosition);
01785         secondary->SetGoodForTrackingFlag();       
01786         secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
01787         secondary->SetWeight(pw[index]);           
01788         fParticleChangeForRadDecay.AddSecondary(secondary); 
01789       }
01790       // make sure the original track is set to stop and its kinematic energy collected
01791       // 
01792       //theTrack.SetTrackStatus(fStopButAlive);
01793       //energyDeposit += theParticle->GetKineticEnergy();
01794 
01795     } // End of Variance Reduction 
01796 
01797     // Kill the parent particle
01798     fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
01799     fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit); 
01800     fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
01801     // Reset NumberOfInteractionLengthLeft.
01802     ClearNumberOfInteractionLengthLeft();
01803 
01804     return &fParticleChangeForRadDecay ;
01805   }
01806 } 
01807 
01808 
01809 G4DecayProducts*
01810 G4RadioactiveDecay::DoDecay(G4ParticleDefinition& theParticleDef)
01811 {
01812   G4DecayProducts* products = 0;
01813 
01814   // follow the decaytable and generate the secondaries...
01815 #ifdef G4VERBOSE
01816   if (GetVerboseLevel()>0) G4cout<<"Begin of DoDecay..."<<G4endl;
01817 #endif
01818 
01819   G4DecayTable* theDecayTable = theParticleDef.GetDecayTable();
01820 
01821   // Choose a decay channel.
01822 #ifdef G4VERBOSE
01823   if (GetVerboseLevel()>0) G4cout <<"Selecte a channel..."<<G4endl;
01824 #endif
01825 
01826   G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel();
01827   if (theDecayChannel == 0) {
01828     // Decay channel not found.
01829     G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
01830     G4cerr <<G4endl;
01831     theDecayTable ->DumpInfo();
01832   } else {
01833     // A decay channel has been identified, so execute the DecayIt.
01834 #ifdef G4VERBOSE
01835     if (GetVerboseLevel()>1) {
01836       G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
01837       G4cerr <<theDecayChannel <<G4endl;
01838     }
01839 #endif
01840     G4double tempmass = theParticleDef.GetPDGMass();
01841     products = theDecayChannel->DecayIt(tempmass);
01842     // Apply directional bias if requested by user
01843     CollimateDecay(products);
01844   }
01845 
01846   return products;
01847 }
01848 
01849 
01850 // Apply directional bias for "visible" daughters (e+-, gamma, n, alpha)
01851 
01852 void G4RadioactiveDecay::CollimateDecay(G4DecayProducts* products) {
01853   if (origin == forceDecayDirection) return;    // No collimation requested
01854   if (180.*deg == forceDecayHalfAngle) return;
01855   if (0 == products || 0 == products->entries()) return;
01856 
01857 #ifdef G4VERBOSE
01858   if (GetVerboseLevel()>0) G4cout<<"Begin of CollimateDecay..."<<G4endl;
01859 #endif
01860 
01861   // Particles suitable for directional biasing (for if-blocks below)
01862   static const G4ParticleDefinition* electron = G4Electron::Definition();
01863   static const G4ParticleDefinition* positron = G4Positron::Definition();
01864   static const G4ParticleDefinition* neutron  = G4Neutron::Definition();
01865   static const G4ParticleDefinition* gamma    = G4Gamma::Definition();
01866   static const G4ParticleDefinition* alpha    = G4Alpha::Definition();
01867 
01868   G4ThreeVector newDirection;           // Re-use to avoid memory churn
01869   for (G4int i=0; i<products->entries(); i++) {
01870     G4DynamicParticle* daughter = (*products)[i];
01871     const G4ParticleDefinition* daughterType = daughter->GetParticleDefinition();
01872 
01873     if (daughterType == electron || daughterType == positron ||
01874         daughterType == neutron || daughterType == gamma ||
01875         daughterType == alpha) CollimateDecayProduct(daughter);
01876   }
01877 }
01878 
01879 void G4RadioactiveDecay::CollimateDecayProduct(G4DynamicParticle* daughter) {
01880 #ifdef G4VERBOSE
01881   if (GetVerboseLevel() > 1) {
01882     G4cout << "CollimateDecayProduct for daughter "
01883            << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
01884   }
01885 #endif
01886 
01887   G4ThreeVector collimate = ChooseCollimationDirection();
01888   if (origin != collimate) daughter->SetMomentumDirection(collimate);
01889 }
01890 
01891 
01892 // Choose random direction within collimation cone
01893 
01894 G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection() const {
01895   if (origin == forceDecayDirection) return origin;     // Don't do collimation
01896   if (forceDecayHalfAngle == 180.*deg) return origin;
01897 
01898   G4ThreeVector dir = forceDecayDirection;
01899 
01900   // Return direction offset by random throw
01901   if (forceDecayHalfAngle > 0.) {
01902     // Generate uniform direction around central axis
01903     G4double phi = 2.*pi*G4UniformRand();
01904     G4double cosMin = std::cos(forceDecayHalfAngle);
01905     G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin;   // [cosMin,1.)
01906     
01907     dir.setPhi(dir.phi()+phi);
01908     dir.setTheta(dir.theta()+std::acos(cosTheta));
01909   }
01910 
01911 #ifdef G4VERBOSE
01912   if (GetVerboseLevel()>1)
01913     G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
01914 #endif
01915 
01916   return dir;
01917 }

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