G4RadioactiveDecay Class Reference

#include <G4RadioactiveDecay.hh>

Inheritance diagram for G4RadioactiveDecay:

Public Member Functions
	G4RadioactiveDecay (const G4String &processName="RadioactiveDecay")
	~G4RadioactiveDecay ()
G4bool	IsApplicable (const G4ParticleDefinition &)
G4bool	IsLoaded (const G4ParticleDefinition &)
void	SelectAVolume (const G4String aVolume)
void	DeselectAVolume (const G4String aVolume)
void	SelectAllVolumes ()
void	DeselectAllVolumes ()
void	SetDecayBias (G4String filename)
void	SetHLThreshold (G4double hl)
void	SetICM (G4bool icm)
void	SetARM (G4bool arm)
void	SetSourceTimeProfile (G4String filename)
G4bool	IsRateTableReady (const G4ParticleDefinition &)
void	AddDecayRateTable (const G4ParticleDefinition &)
void	GetDecayRateTable (const G4ParticleDefinition &)
void	SetDecayRate (G4int, G4int, G4double, G4int, std::vector< G4double >, std::vector< G4double >)
std::vector< G4RadioactivityTable * >	GetTheRadioactivityTables ()
G4DecayTable *	LoadDecayTable (G4ParticleDefinition &theParentNucleus)
void	AddUserDecayDataFile (G4int Z, G4int A, G4String filename)
void	SetVerboseLevel (G4int value)
G4int	GetVerboseLevel () const
void	SetNucleusLimits (G4NucleusLimits theNucleusLimits1)
G4NucleusLimits	GetNucleusLimits () const
void	SetAnalogueMonteCarlo (G4bool r)
void	SetFBeta (G4bool r)
G4bool	IsAnalogueMonteCarlo ()
void	SetBRBias (G4bool r)
void	SetSplitNuclei (G4int r)
G4int	GetSplitNuclei ()
void	SetDecayDirection (const G4ThreeVector &theDir)
const G4ThreeVector &	GetDecayDirection () const
void	SetDecayHalfAngle (G4double halfAngle=0.*CLHEP::deg)
G4double	GetDecayHalfAngle () const
void	SetDecayCollimation (const G4ThreeVector &theDir, G4double halfAngle=0.*CLHEP::deg)
void	BuildPhysicsTable (const G4ParticleDefinition &)
Protected Member Functions
G4VParticleChange *	DecayIt (const G4Track &theTrack, const G4Step &theStep)
G4DecayProducts *	DoDecay (G4ParticleDefinition &theParticleDef)
void	CollimateDecay (G4DecayProducts *products)
void	CollimateDecayProduct (G4DynamicParticle *product)
G4ThreeVector	ChooseCollimationDirection () const
G4double	GetMeanFreePath (const G4Track &theTrack, G4double previousStepSize, G4ForceCondition *condition)
G4double	GetMeanLifeTime (const G4Track &theTrack, G4ForceCondition *condition)
G4double	GetTaoTime (const G4double, const G4double)
G4double	GetDecayTime ()
G4int	GetDecayTimeBin (const G4double aDecayTime)

---

Detailed Description

Definition at line 80 of file G4RadioactiveDecay.hh.

---

Constructor & Destructor Documentation

G4RadioactiveDecay::G4RadioactiveDecay

(

const G4String &

processName = "RadioactiveDecay"

)

Definition at line 144 of file G4RadioactiveDecay.cc.

References fRadioactiveDecay, G4cout, G4endl, G4ParticleTable::GetParticleTable(), GetVerboseLevel(), G4VProcess::pParticleChange, G4IonTable::RegisterIsotopeTable(), SelectAllVolumes(), and G4VProcess::SetProcessSubType().

 : G4VRestDiscreteProcess(processName, fDecay), HighestValue(20.0),
   isInitialised(false), forceDecayDirection(0.,0.,0.), 
   forceDecayHalfAngle(0.*deg), verboseLevel(0)
{
#ifdef G4VERBOSE
  if (GetVerboseLevel()>1) {
    G4cout <<"G4RadioactiveDecay constructor    Name: ";
    G4cout <<processName << G4endl;   }
#endif

  SetProcessSubType(fRadioactiveDecay);

  theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
  theIsotopeTable              = new G4RIsotopeTable();
  pParticleChange              = &fParticleChangeForRadDecay;

  // Now register the Isotope table with G4IonTable.
  G4IonTable *theIonTable =
    (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
  G4VIsotopeTable *aVirtualTable = theIsotopeTable;
  theIonTable->RegisterIsotopeTable(aVirtualTable);

  // Reset the contents of the list of nuclei for which decay scheme data
  // have been loaded.
  LoadedNuclei.clear();

  //
  //Reset the list of user define data file
  //
  theUserRadioactiveDataFiles.clear();

  //
  //
  // Apply default values.
  //
  NSourceBin  = 1;
  SBin[0]     = 0.* s;
  SBin[1]     = 1.* s;
  SProfile[0] = 1.;
  SProfile[1] = 0.;
  NDecayBin   = 1;
  DBin[0]     = 0. * s ;
  DBin[1]     = 1. * s;
  DProfile[0] = 1.;
  DProfile[1] = 0.;
  decayWindows[0] = 0;
  G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
  theRadioactivityTables.push_back(rTable);
  NSplit      = 1;
  AnalogueMC  = true ;
  FBeta       = false ;
  BRBias      = true ;
  applyICM    = true ;
  applyARM    = true ;
  halflifethreshold = -1.*second;
  //
  // RDM applies to xall logical volumes as default
  isAllVolumesMode=true;
  SelectAllVolumes();
}

G4RadioactiveDecay::~G4RadioactiveDecay

(

)

Definition at line 207 of file G4RadioactiveDecay.cc.

{
  delete theRadioactiveDecaymessenger;
}

---

Member Function Documentation

void G4RadioactiveDecay::AddDecayRateTable

(

const G4ParticleDefinition &

)

Definition at line 1064 of file G4RadioactiveDecay.cc.

References G4DecayTable::entries(), G4cout, G4endl, G4ParticleDefinition::GetAtomicMass(), G4VDecayChannel::GetBR(), G4NuclearDecayChannel::GetDaughterExcitation(), G4NuclearDecayChannel::GetDaughterNucleus(), G4DecayTable::GetDecayChannel(), G4NuclearDecayChannel::GetDecayMode(), G4ParticleDefinition::GetDecayTable(), G4NuclearLevelStore::GetInstance(), G4IonTable::GetIon(), G4NuclearLevelStore::GetManager(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGLifeTime(), GetVerboseLevel(), G4DecayTable::Insert(), IsApplicable(), IsLoaded(), IT, LoadDecayTable(), G4NuclearLevelManager::NearestLevel(), ns, G4NuclearLevelManager::NumberOfLevels(), SetDecayRate(), G4ParticleDefinition::SetDecayTable(), G4RadioactiveDecayRateVector::SetIonName(), and G4RadioactiveDecayRateVector::SetItsRates().

Referenced by DecayIt().

{
  // 1) To calculate all the coefficiecies required to derive the
  //    radioactivities for all progeny of theParentNucleus
  //
  // 2) Add the coefficiencies to the decay rate table vector 
  //

  //
  // Create and initialise variables used in the method.
  //
  theDecayRateVector.clear();

  G4int nGeneration = 0;
  std::vector<G4double> rates;
  std::vector<G4double> taos;

  // start rate is -1.
  // Eq.4.26 of the Technical Note
  rates.push_back(-1.);
  //
  //
  G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
  G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
  G4double tao = theParentNucleus.GetPDGLifeTime();
  if (tao < 0.) tao = 1e-100;
  taos.push_back(tao);
  G4int nEntry = 0;

  //fill the decay rate with the intial isotope data
  SetDecayRate(Z,A,E,nGeneration,rates,taos);

  // store the decay rate in decay rate vector
  theDecayRateVector.push_back(theDecayRate);
  nEntry++;

  // now start treating the sencondary generations..

  G4bool stable = false;
  G4int i;
  G4int j;
  G4VDecayChannel* theChannel = 0;
  G4NuclearDecayChannel* theNuclearDecayChannel = 0;
  G4ITDecayChannel* theITChannel = 0;
  G4BetaMinusDecayChannel *theBetaMinusChannel = 0;
  G4BetaPlusDecayChannel *theBetaPlusChannel = 0;
  G4AlphaDecayChannel *theAlphaChannel = 0;
  G4RadioactiveDecayMode theDecayMode;
  G4double theBR = 0.0;
  G4int AP = 0;
  G4int ZP = 0;
  G4int AD = 0;
  G4int ZD = 0;
  G4double EP = 0.;
  std::vector<G4double> TP;
  std::vector<G4double> RP;
  G4ParticleDefinition *theDaughterNucleus;
  G4double daughterExcitation;
  G4ParticleDefinition *aParentNucleus;
  G4IonTable* theIonTable;
  G4DecayTable *aTempDecayTable;
  G4double theRate;
  G4double TaoPlus;
  G4int nS = 0;
  G4int nT = nEntry;
  G4double brs[7];
  //
  theIonTable =
    (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());

  while (!stable) {
    nGeneration++;
    for (j = nS; j< nT; j++) {
      ZP = theDecayRateVector[j].GetZ();
      AP = theDecayRateVector[j].GetA();
      EP = theDecayRateVector[j].GetE();
      RP = theDecayRateVector[j].GetDecayRateC();
      TP = theDecayRateVector[j].GetTaos();      
      if (GetVerboseLevel()>0){
        G4cout <<"G4RadioactiveDecay::AddDecayRateTable : "
               << " daughters of ("<< ZP <<", "<<AP<<", "
               << EP <<") "
               << " are being calculated. "       
               <<" generation = "
               << nGeneration << G4endl;
      }
      aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
      if (!IsLoaded(*aParentNucleus)){
        aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
      }
        
      G4DecayTable* theDecayTable = new G4DecayTable();
      aTempDecayTable = aParentNucleus->GetDecayTable();
      for (i=0; i< 7; i++) brs[i] = 0.0;

      //
      // Go through the decay table and to combine the same decay channels
      //
      for (i=0; i<aTempDecayTable->entries(); i++) {
        theChannel             = aTempDecayTable->GetDecayChannel(i);
        theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
        theDecayMode           = theNuclearDecayChannel->GetDecayMode();
        daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation ();
        theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus () ;
        AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
        ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();  
        G4NuclearLevelManager* levelManager = 
           G4NuclearLevelStore::GetInstance()->GetManager(ZD, AD);
        if (levelManager->NumberOfLevels() ) {
          const G4NuclearLevel* level =
              levelManager->NearestLevel (daughterExcitation);

          if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
            // Level half-life is in ns and the threshold is set to 1 micros by default, user can set it via the UI command
            if (level->HalfLife()*ns >= halflifethreshold ){    
              // save the metastable nucleus 
              theDecayTable->Insert(theChannel);
            } 
            else{
              brs[theDecayMode] += theChannel->GetBR();
            }
          }
          else {
            brs[theDecayMode] += theChannel->GetBR();
          }
        }
        else{
          brs[theDecayMode] += theChannel->GetBR();
        }
      }     
      brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
      brs[3] = brs[4] =brs[5] =  0.0;
      for (i= 0; i<7; i++){
        if (brs[i] > 0.) {
          switch ( i ) {
          case 0:
            //
            //
            // Decay mode is isomeric transition.
            //

            theITChannel =  new G4ITDecayChannel
              (0, (const G4Ions*) aParentNucleus, brs[0]);

            theDecayTable->Insert(theITChannel);
            break;

          case 1:
            //
            //
            // Decay mode is beta-.
            //
            theBetaMinusChannel = new G4BetaMinusDecayChannel (0, aParentNucleus,
                                                               brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
            theDecayTable->Insert(theBetaMinusChannel);

            break;

          case 2:
            //
            //
            // Decay mode is beta+ + EC.
            //
            theBetaPlusChannel = new G4BetaPlusDecayChannel (GetVerboseLevel(), aParentNucleus,
                                                             brs[2], 0.*MeV, 0.*MeV, 1, false, 0);
            theDecayTable->Insert(theBetaPlusChannel);
            break;                    

          case 6:
            //
            //
            // Decay mode is alpha.
            //
            theAlphaChannel = new G4AlphaDecayChannel(GetVerboseLevel(),
                                                      aParentNucleus,
                                                      brs[6], 0.*MeV, 0.*MeV);
            theDecayTable->Insert(theAlphaChannel);
            break;

          default:
            break;
          }
        }
      }
      // 
      // loop over all branches in theDecayTable
      //
      for (i = 0; i < theDecayTable->entries(); i++){
        theChannel = theDecayTable->GetDecayChannel(i);
        theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
        theBR = theChannel->GetBR();
        theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
        //  first check if the decay of the original nucleus is an IT channel, if true create a new groud-level nucleus
        if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
          A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
          Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
          theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
        }
        if (IsApplicable(*theDaughterNucleus) &&
            theBR && 
            aParentNucleus != theDaughterNucleus) { 
          // need to make sure daugher has decaytable
          if (!IsLoaded(*theDaughterNucleus)){
            theDaughterNucleus->SetDecayTable(LoadDecayTable(*theDaughterNucleus));
          }
          if (theDaughterNucleus->GetDecayTable()->entries() ) {
            //
            A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
            Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
            E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();

            TaoPlus = theDaughterNucleus->GetPDGLifeTime();
            //          cout << TaoPlus <<G4endl;

            if (TaoPlus <= 0.)  TaoPlus = 1e-100;



            // first set the taos, one simply need to add to the parent ones
            taos.clear();
            taos = TP;
            size_t k;
            //check that TaoPlus differs from other taos from at least 1.e5 relative difference
            //for (k = 0; k < TP.size(); k++){
            //  if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
            //}
            taos.push_back(TaoPlus);
            // now calculate the coefficiencies
            //
            // they are in two parts, first the less than n ones
            // Eq 4.24 of the TN
            rates.clear();
            long double ta1,ta2;
            ta2 = (long double)TaoPlus;
            for (k = 0; k < RP.size(); k++){
              ta1 = (long double)TP[k];
              if (ta1 == ta2) {
                theRate = 1.e100;
              }else{
                theRate = ta1/(ta1-ta2);}
              theRate = theRate * theBR * RP[k];
              rates.push_back(theRate);
            }
            //
            // the sencond part: the n:n coefficiency
            // Eq 4.25 of the TN.  Note Yn+1 is zero apart from Y1 which is -1 as treated at line 1013
            // 
            theRate = 0.;
            long double aRate, aRate1;
            aRate1 = 0.L;
            for (k = 0; k < RP.size(); k++){
              ta1 = (long double)TP[k];
              if (ta1 == ta2 ) {
                aRate = 1.e100;
              }else {
                aRate = ta2/(ta1-ta2);}
              aRate = aRate * (long double)(theBR * RP[k]);
              aRate1 += aRate;
            }
            theRate = -aRate1;
            rates.push_back(theRate);         
            SetDecayRate (Z,A,E,nGeneration,rates,taos);
            theDecayRateVector.push_back(theDecayRate);
            nEntry++;
          }
        } // end of testing daughter nucleus
      } // end of i loop( the branches) 
      //      delete theDecayTable;

    } //end of for j loop
    nS = nT;
    nT = nEntry;
    if (nS == nT) stable = true;
  }

  //end of while loop
  // the calculation completed here


  // fill the first part of the decay rate table
  // which is the name of the original particle (isotope) 
  //
  theDecayRateTable.SetIonName(theParentNucleus.GetParticleName()); 
  //
  //
  // now fill the decay table with the newly completed decay rate vector
  //

  theDecayRateTable.SetItsRates(theDecayRateVector);

  //
  // finally add the decayratetable to the tablevector
  //
  theDecayRateTableVector.push_back(theDecayRateTable);
}

void G4RadioactiveDecay::AddUserDecayDataFile	(	G4int	Z,
		G4int	A,
		G4String	filename
	)

Definition at line 1028 of file G4RadioactiveDecay.cc.

References G4RIsotopeTable::AddUserDecayDataFile(), G4cout, and G4endl.

Referenced by G4RadioactiveDecaymessenger::SetNewValue().

{ if (Z<1 || A<2) {
        G4cout<<"Z and A not valid!"<<G4endl;
  }

  std::ifstream DecaySchemeFile(filename);
  if (DecaySchemeFile){
        G4int ID_ion=A*1000+Z;
        theUserRadioactiveDataFiles[ID_ion]=filename;
        theIsotopeTable->AddUserDecayDataFile(Z,A,filename);
  }
  else {
        G4cout<<"The file "<<filename<<" does not exist!"<<G4endl;
  }
}

void G4RadioactiveDecay::BuildPhysicsTable

(

const G4ParticleDefinition &

)

[virtual]

Reimplemented from G4VProcess.

Definition at line 670 of file G4RadioactiveDecay.cc.

References G4LossTableManager::AtomDeexcitation(), G4VAtomDeexcitation::InitialiseAtomicDeexcitation(), and G4LossTableManager::Instance().

{
  if(!isInitialised) {
    isInitialised = true;
    G4VAtomDeexcitation* p = G4LossTableManager::Instance()->AtomDeexcitation();
    if(p) { p->InitialiseAtomicDeexcitation(); }
  }
}

G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection

(

)

const [protected]

Definition at line 1894 of file G4RadioactiveDecay.cc.

References G4cout, G4endl, G4UniformRand, GetVerboseLevel(), and G4INCL::Math::pi.

Referenced by CollimateDecayProduct().

                                                                   {
  if (origin == forceDecayDirection) return origin;     // Don't do collimation
  if (forceDecayHalfAngle == 180.*deg) return origin;

  G4ThreeVector dir = forceDecayDirection;

  // Return direction offset by random throw
  if (forceDecayHalfAngle > 0.) {
    // Generate uniform direction around central axis
    G4double phi = 2.*pi*G4UniformRand();
    G4double cosMin = std::cos(forceDecayHalfAngle);
    G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin;   // [cosMin,1.)
    
    dir.setPhi(dir.phi()+phi);
    dir.setTheta(dir.theta()+std::acos(cosTheta));
  }

#ifdef G4VERBOSE
  if (GetVerboseLevel()>1)
    G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
#endif

  return dir;
}

void G4RadioactiveDecay::CollimateDecay

(

G4DecayProducts *

products

)

[protected]

Definition at line 1852 of file G4RadioactiveDecay.cc.

References G4InuclParticleNames::alpha, CollimateDecayProduct(), G4Alpha::Definition(), G4Gamma::Definition(), G4Neutron::Definition(), G4Positron::Definition(), G4Electron::Definition(), G4DecayProducts::entries(), G4cout, G4endl, GetVerboseLevel(), and neutron.

Referenced by DoDecay().

                                                                 {
  if (origin == forceDecayDirection) return;    // No collimation requested
  if (180.*deg == forceDecayHalfAngle) return;
  if (0 == products || 0 == products->entries()) return;

#ifdef G4VERBOSE
  if (GetVerboseLevel()>0) G4cout<<"Begin of CollimateDecay..."<<G4endl;
#endif

  // Particles suitable for directional biasing (for if-blocks below)
  static const G4ParticleDefinition* electron = G4Electron::Definition();
  static const G4ParticleDefinition* positron = G4Positron::Definition();
  static const G4ParticleDefinition* neutron  = G4Neutron::Definition();
  static const G4ParticleDefinition* gamma    = G4Gamma::Definition();
  static const G4ParticleDefinition* alpha    = G4Alpha::Definition();

  G4ThreeVector newDirection;           // Re-use to avoid memory churn
  for (G4int i=0; i<products->entries(); i++) {
    G4DynamicParticle* daughter = (*products)[i];
    const G4ParticleDefinition* daughterType = daughter->GetParticleDefinition();

    if (daughterType == electron || daughterType == positron ||
        daughterType == neutron || daughterType == gamma ||
        daughterType == alpha) CollimateDecayProduct(daughter);
  }
}

void G4RadioactiveDecay::CollimateDecayProduct

(

G4DynamicParticle *

product

)

[protected]

Definition at line 1879 of file G4RadioactiveDecay.cc.

References ChooseCollimationDirection(), G4cout, G4endl, G4DynamicParticle::GetParticleDefinition(), G4ParticleDefinition::GetParticleName(), GetVerboseLevel(), and G4DynamicParticle::SetMomentumDirection().

Referenced by CollimateDecay().

                                                                          {
#ifdef G4VERBOSE
  if (GetVerboseLevel() > 1) {
    G4cout << "CollimateDecayProduct for daughter "
           << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
  }
#endif

  G4ThreeVector collimate = ChooseCollimationDirection();
  if (origin != collimate) daughter->SetMomentumDirection(collimate);
}

G4VParticleChange * G4RadioactiveDecay::DecayIt	(	const G4Track &	theTrack,
		const G4Step &	theStep
	)			[protected]

Definition at line 1457 of file G4RadioactiveDecay.cc.

References AddDecayRateTable(), G4ParticleChangeForRadDecay::AddSecondary(), G4DecayProducts::Boost(), G4VProcess::ClearNumberOfInteractionLengthLeft(), DoDecay(), G4DecayTable::DumpInfo(), G4DecayProducts::DumpInfo(), G4DecayProducts::entries(), G4DecayTable::entries(), fStopAndKill, fStopButAlive, G4cerr, G4cout, G4endl, G4UniformRand, G4ParticleDefinition::GetBaryonNumber(), G4DecayTable::GetDecayChannel(), GetDecayRateTable(), G4ParticleDefinition::GetDecayTable(), GetDecayTime(), GetDecayTimeBin(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4Track::GetGlobalTime(), G4IonTable::GetIon(), G4DynamicParticle::GetKineticEnergy(), G4Track::GetLocalTime(), G4VPhysicalVolume::GetLogicalVolume(), G4DynamicParticle::GetMomentumDirection(), G4LogicalVolume::GetName(), G4DynamicParticle::GetParticleDefinition(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGMass(), G4Track::GetPosition(), GetTaoTime(), G4Track::GetTouchableHandle(), G4Track::GetTrackStatus(), GetVerboseLevel(), G4Track::GetVolume(), G4Track::GetWeight(), G4ParticleChangeForDecay::Initialize(), IsApplicable(), G4DecayProducts::IsChecked(), IsLoaded(), IsRateTableReady(), LoadDecayTable(), CLHEP::detail::n, ns, G4DecayProducts::PopProducts(), G4VParticleChange::ProposeLocalEnergyDeposit(), G4ParticleChangeForDecay::ProposeLocalTime(), G4VParticleChange::ProposeTrackStatus(), G4ParticleDefinition::SetDecayTable(), and G4VParticleChange::SetNumberOfSecondaries().

{
  // Initialize the G4ParticleChange object. Get the particle details and the
  // decay table.



  fParticleChangeForRadDecay.Initialize(theTrack);
  const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();


  G4ParticleDefinition *theParticleDef = theParticle->GetDefinition();


  // First check whether RDM applies to the current logical volume
  if (!isAllVolumesMode){
   if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
                          theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
#ifdef G4VERBOSE
    if (GetVerboseLevel()>0) {
      G4cout <<"G4RadioactiveDecay::DecayIt : "
             << theTrack.GetVolume()->GetLogicalVolume()->GetName()
             << " is not selected for the RDM"<< G4endl;
      G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
      G4cout << " The Valid volumes are " << G4endl;
      for (size_t i = 0; i< ValidVolumes.size(); i++) G4cout << ValidVolumes[i] << G4endl;
    }
#endif
    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);

    // Kill the parent particle.

    fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
    ClearNumberOfInteractionLengthLeft();
    return &fParticleChangeForRadDecay;
   }
  }

  // now check is the particle is valid for RDM

  if (!(IsApplicable(*theParticleDef))) { 
    //
    // The particle is not a Ion or outside the nucleuslimits for decay
    //
#ifdef G4VERBOSE
    if (GetVerboseLevel()>0) {
      G4cerr <<"G4RadioactiveDecay::DecayIt : "
             <<theParticleDef->GetParticleName() 
             << " is not a valid nucleus for the RDM"<< G4endl;
    }
#endif
    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);

    //
    // Kill the parent particle.
    //
    fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
    ClearNumberOfInteractionLengthLeft();
    return &fParticleChangeForRadDecay;
  }

  if (!IsLoaded(*theParticleDef))
    theParticleDef->SetDecayTable(LoadDecayTable(*theParticleDef));

  G4DecayTable* theDecayTable = theParticleDef->GetDecayTable();

  if (theDecayTable == 0 || theDecayTable->entries() == 0) {
      // There are no data in the decay table.  Set the particle change parameters
      // to indicate this.
#ifdef G4VERBOSE
      if (GetVerboseLevel()>0)
        {
          G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined  for ";
          G4cerr <<theParticleDef->GetParticleName() <<G4endl;
        }
#endif
      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);

      // Kill the parent particle.
      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
      ClearNumberOfInteractionLengthLeft();
      return &fParticleChangeForRadDecay;

  } else { 
    // Data found.  Try to decay nucleus
    G4double energyDeposit = 0.0;
    G4double finalGlobalTime = theTrack.GetGlobalTime();
    G4double finalLocalTime = theTrack.GetLocalTime();
    G4int index;
    G4ThreeVector currentPosition;
    currentPosition = theTrack.GetPosition();

    // Check whether use Analogue or VR implementation
    if (AnalogueMC) {
#ifdef G4VERBOSE
      if (GetVerboseLevel() > 0)
        G4cout <<"DecayIt:  Analogue MC version " << G4endl;
#endif

      G4DecayProducts* products = DoDecay(*theParticleDef);

      // Check if the product is the same as input and kill the track if
      // necessary to prevent infinite loop (11/05/10, F.Lei)
      if ( products->entries() == 1) {
        fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
        fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
        fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
        ClearNumberOfInteractionLengthLeft();
        return &fParticleChangeForRadDecay;
      }

      // Get parent particle information and boost the decay products to the
      // laboratory frame based on this information.


      //The Parent Energy used for the boost should be the total energy of
      // the nucleus of the parent ion without the energy of the shell electrons
      // (correction for bug 1359 by L. Desorgher)
      G4double ParentEnergy =  theParticle->GetKineticEnergy()+theParticle->GetParticleDefinition()->GetPDGMass();
      G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());



      if (theTrack.GetTrackStatus() == fStopButAlive) { //this condition seems to be always True, further investigation is needed (L.Desorgher)

        // The particle is decayed at rest.
        // since the time is still for rest particle in G4 we need to add the
        // additional time lapsed between the particle come to rest and the
        // actual decay.  This time is simply sampled with the mean-life of
        // the particle.  But we need to protect the case PDGTime < 0.
        // (F.Lei 11/05/10)
        G4double temptime = -std::log( G4UniformRand())
                            *theParticleDef->GetPDGLifeTime();
        if (temptime < 0.) temptime = 0.; 
        finalGlobalTime += temptime;
        finalLocalTime += temptime;
        energyDeposit += theParticle->GetKineticEnergy();
      }
      products->Boost( ParentEnergy, ParentDirection);



      // Add products in theParticleChangeForRadDecay.
      G4int numberOfSecondaries = products->entries();
      fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
#ifdef G4VERBOSE
      if (GetVerboseLevel()>1) {
        G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
        G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
        G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
        G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
        G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
        G4cout << G4endl;
        G4cout <<"G4Decay::DecayIt  : decay products in Lab. Frame" <<G4endl;
        products->DumpInfo();
        products->IsChecked();
      }
#endif
      for (index=0; index < numberOfSecondaries; index++) {
        G4Track* secondary = new G4Track(products->PopProducts(),
                                         finalGlobalTime, currentPosition);
        secondary->SetGoodForTrackingFlag();
        secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
        fParticleChangeForRadDecay.AddSecondary(secondary);
      }
      delete products;
      // end of analogue MC algorithm

    } else {
      // Variance Reduction Method
#ifdef G4VERBOSE
      if (GetVerboseLevel()>0)
        G4cout << "DecayIt: Variance Reduction version " << G4endl;
#endif
      if (!IsRateTableReady(*theParticleDef)) {
        // if the decayrates are not ready, calculate them and 
        // add to the rate table vector 
        AddDecayRateTable(*theParticleDef);
      }
      //retrieve the rates 
      GetDecayRateTable(*theParticleDef);

      // declare some of the variables required in the implementation
      G4ParticleDefinition* parentNucleus;
      G4IonTable* theIonTable;
      G4int PZ;
      G4int PA;
      G4double PE;
      std::vector<G4double> PT;
      std::vector<G4double> PR;
      G4double taotime;
      long double decayRate;

      size_t i;
      size_t j;
      G4int numberOfSecondaries;
      G4int totalNumberOfSecondaries = 0;
      G4double currentTime = 0.;
      G4int ndecaych;
      G4DynamicParticle* asecondaryparticle;
      std::vector<G4DynamicParticle*> secondaryparticles;
      std::vector<G4double> pw;
      std::vector<G4double> ptime;
      pw.clear();
      ptime.clear();

      //now apply the nucleus splitting
      for (G4int n = 0; n < NSplit; n++) {
        // Get the decay time following the decay probability function 
        // suppllied by user  
        G4double theDecayTime = GetDecayTime();
        G4int nbin = GetDecayTimeBin(theDecayTime);
            
        // calculate the first part of the weight function  
        G4double weight1 = 1.; 
        if (nbin == 1) {
          weight1 = 1./DProfile[nbin-1] 
                    *(DBin[nbin]-DBin[nbin-1])/NSplit;
        } else if (nbin > 1) {
          weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
                    *(DBin[nbin]-DBin[nbin-1])/NSplit;
        }

        // it should be calculated in seconds
        weight1 /= s ;
            
        // loop over all the possible secondaries of the nucleus
        // the first one is itself.
        for (i = 0; i < theDecayRateVector.size(); i++) {
          PZ = theDecayRateVector[i].GetZ();
          PA = theDecayRateVector[i].GetA();
          PE = theDecayRateVector[i].GetE();
          PT = theDecayRateVector[i].GetTaos();
          PR = theDecayRateVector[i].GetDecayRateC();

          // Calculate the decay rate of the isotope
          // decayRate is the radioactivity of isotope (PZ,PA,PE) at the 
          // time 'theDecayTime'
          // it will be used to calculate the statistical weight of the 
          // decay products of this isotope

          //  G4cout <<"PA= "<< PA << " PZ= " << PZ << " PE= "<< PE  <<G4endl;
          decayRate = 0.L;
          for (j = 0; j < PT.size(); j++) {
            taotime = GetTaoTime(theDecayTime,PT[j]);
            decayRate -= PR[j] * (long double)taotime;
            // Eq.4.23 of of the TN
            // note the negative here is required as the rate in the
            // equation is defined to be negative, 
            // i.e. decay away, but we need positive value here.

            // G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t"
            //        << decayRate << G4endl;           
          }

          // add the isotope to the radioactivity tables
          //  G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
          //  G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
          theRadioactivityTables[decayWindows[nbin-1]]->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());

          // Now calculate the statistical weight
          // One needs to fold the source bias function with the decaytime
          // also need to include the track weight! (F.Lei, 28/10/10)
          G4double weight = weight1*decayRate*theTrack.GetWeight();



          // decay the isotope 
          theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
          parentNucleus = theIonTable->GetIon(PZ,PA,PE);

          // Create a temprary products buffer.
          // Its contents to be transfered to the products at the end of the loop
          G4DecayProducts* tempprods = 0;

          // Decide whether to apply branching ratio bias or not             
          if (BRBias) {
            G4DecayTable* decayTable = parentNucleus->GetDecayTable();
            ndecaych = G4int(decayTable->entries()*G4UniformRand());
            G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
            if (theDecayChannel == 0) {
              // Decay channel not found.
#ifdef G4VERBOSE
              if (GetVerboseLevel()>0) {
                G4cerr << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
                G4cerr << " for this nucleus; decay as if no biasing active ";
                G4cerr << G4endl;
                decayTable ->DumpInfo();
              }
#endif
              tempprods = DoDecay(*parentNucleus);  // DHW 6 Dec 2010 - do decay as if no biasing
                                                    //           to avoid deref of temppprods = 0
            } else {
              // A decay channel has been identified, so execute the DecayIt.
              G4double tempmass = parentNucleus->GetPDGMass();
              tempprods = theDecayChannel->DecayIt(tempmass);
              weight *= (theDecayChannel->GetBR())*(decayTable->entries());
            }
          } else {
            tempprods = DoDecay(*parentNucleus);
          }

          // save the secondaries for buffers
          numberOfSecondaries = tempprods->entries();
          currentTime = finalGlobalTime + theDecayTime;
          for (index = 0; index < numberOfSecondaries; index++) {
            asecondaryparticle = tempprods->PopProducts();
            if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5) {
              pw.push_back(weight);
              ptime.push_back(currentTime);
              secondaryparticles.push_back(asecondaryparticle);
            }
          }
          delete tempprods;

        } // end of i loop
      } // end of n loop 

      // now deal with the secondaries in the two stl containers
      // and submmit them back to the tracking manager
      totalNumberOfSecondaries = pw.size();
      fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
      for (index=0; index < totalNumberOfSecondaries; index++) { 
        G4Track* secondary = new G4Track(secondaryparticles[index],
                                         ptime[index], currentPosition);
        secondary->SetGoodForTrackingFlag();       
        secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
        secondary->SetWeight(pw[index]);           
        fParticleChangeForRadDecay.AddSecondary(secondary); 
      }
      // make sure the original track is set to stop and its kinematic energy collected
      // 
      //theTrack.SetTrackStatus(fStopButAlive);
      //energyDeposit += theParticle->GetKineticEnergy();

    } // End of Variance Reduction 

    // Kill the parent particle
    fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit); 
    fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
    // Reset NumberOfInteractionLengthLeft.
    ClearNumberOfInteractionLengthLeft();

    return &fParticleChangeForRadDecay ;
  }
} 

void G4RadioactiveDecay::DeselectAllVolumes

(

)

Definition at line 324 of file G4RadioactiveDecay.cc.

References G4cout, G4endl, and GetVerboseLevel().

{
  ValidVolumes.clear();
  isAllVolumesMode=false;
#ifdef G4VERBOSE
  if (GetVerboseLevel()>0)
    G4cout << " RDM removed from all volumes" << G4endl; 
#endif
}

void G4RadioactiveDecay::DeselectAVolume

(

const G4String

aVolume

)

Definition at line 271 of file G4RadioactiveDecay.cc.

References G4cerr, G4cout, G4endl, G4LogicalVolumeStore::GetInstance(), G4LogicalVolume::GetName(), and GetVerboseLevel().

{
  G4LogicalVolumeStore *theLogicalVolumes;
  G4LogicalVolume *volume;
  theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
    volume=(*theLogicalVolumes)[i];
    if (volume->GetName() == aVolume) {
      std::vector<G4String>::iterator location;
      location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
      if (location != ValidVolumes.end()) {
        ValidVolumes.erase(location);
        std::sort(ValidVolumes.begin(), ValidVolumes.end());
        isAllVolumesMode =false;
      } else {
        G4cerr << " DeselectVolume:" << aVolume << " is not in the list"<< G4endl; 
      }   
#ifdef G4VERBOSE
      if (GetVerboseLevel()>0)
        G4cout << " DeselectVolume: " << aVolume << " is removed from list"<<G4endl; 
#endif
    } else if (i ==  theLogicalVolumes->size()) {
      G4cerr << " DeselectVolume:" << aVolume
             << "is not a valid logical volume name"<< G4endl; 
    }
  }
}

G4DecayProducts * G4RadioactiveDecay::DoDecay

(

G4ParticleDefinition &

theParticleDef

)

[protected]

Definition at line 1810 of file G4RadioactiveDecay.cc.

References CollimateDecay(), G4VDecayChannel::DecayIt(), G4DecayTable::DumpInfo(), G4cerr, G4cout, G4endl, G4ParticleDefinition::GetDecayTable(), G4ParticleDefinition::GetPDGMass(), GetVerboseLevel(), and G4DecayTable::SelectADecayChannel().

Referenced by DecayIt().

{
  G4DecayProducts* products = 0;

  // follow the decaytable and generate the secondaries...
#ifdef G4VERBOSE
  if (GetVerboseLevel()>0) G4cout<<"Begin of DoDecay..."<<G4endl;
#endif

  G4DecayTable* theDecayTable = theParticleDef.GetDecayTable();

  // Choose a decay channel.
#ifdef G4VERBOSE
  if (GetVerboseLevel()>0) G4cout <<"Selecte a channel..."<<G4endl;
#endif

  G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel();
  if (theDecayChannel == 0) {
    // Decay channel not found.
    G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
    G4cerr <<G4endl;
    theDecayTable ->DumpInfo();
  } else {
    // A decay channel has been identified, so execute the DecayIt.
#ifdef G4VERBOSE
    if (GetVerboseLevel()>1) {
      G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
      G4cerr <<theDecayChannel <<G4endl;
    }
#endif
    G4double tempmass = theParticleDef.GetPDGMass();
    products = theDecayChannel->DecayIt(tempmass);
    // Apply directional bias if requested by user
    CollimateDecay(products);
  }

  return products;
}

const G4ThreeVector& G4RadioactiveDecay::GetDecayDirection

(

)

const [inline]

Definition at line 217 of file G4RadioactiveDecay.hh.

                                                          {
      return forceDecayDirection; 
    }

G4double G4RadioactiveDecay::GetDecayHalfAngle

(

)

const [inline]

Definition at line 225 of file G4RadioactiveDecay.hh.

{return forceDecayHalfAngle;}

void G4RadioactiveDecay::GetDecayRateTable

(

const G4ParticleDefinition &

)

Definition at line 352 of file G4RadioactiveDecay.cc.

References G4cout, G4endl, G4ParticleDefinition::GetParticleName(), and GetVerboseLevel().

Referenced by DecayIt().

{
  G4String aParticleName = aParticle.GetParticleName();

  for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
    if (theDecayRateTableVector[i].GetIonName() == aParticleName) {
      theDecayRateVector = theDecayRateTableVector[i].GetItsRates();
    }
  }
#ifdef G4VERBOSE
  if (GetVerboseLevel() > 0) {
    G4cout << "The DecayRate Table for "
           << aParticleName << " is selected." <<  G4endl;
  }
#endif
}

G4double G4RadioactiveDecay::GetDecayTime

(

)

[protected]

Definition at line 522 of file G4RadioactiveDecay.cc.

References G4cout, G4endl, G4UniformRand, and GetVerboseLevel().

Referenced by DecayIt().

{
  G4double decaytime = 0.;
  G4double rand = G4UniformRand();
  G4int i = 0;
  while ( DProfile[i] < rand) i++;
  rand = G4UniformRand();
  decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
#ifdef G4VERBOSE
  if (GetVerboseLevel()>1)
    {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;}
#endif
  return  decaytime;        
}

G4int G4RadioactiveDecay::GetDecayTimeBin

(

const G4double

aDecayTime

)

[protected]

Definition at line 538 of file G4RadioactiveDecay.cc.

Referenced by DecayIt().

{
  G4int i = 0;
  while ( aDecayTime > DBin[i] ) i++;
  return  i;
}

G4double G4RadioactiveDecay::GetMeanFreePath	(	const G4Track &	theTrack,
		G4double	previousStepSize,
		G4ForceCondition *	condition
	)			[protected, virtual]

Implements G4VRestDiscreteProcess.

Definition at line 593 of file G4RadioactiveDecay.cc.

References DBL_MAX, DBL_MIN, G4cerr, G4cout, G4endl, G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGLifeTime(), G4ParticleDefinition::GetPDGStable(), G4DynamicParticle::GetTotalMomentum(), and GetVerboseLevel().

{
  // get particle
  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();

  // returns the mean free path in GEANT4 internal units
  G4double pathlength;
  G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
  G4double aCtau = c_light * aParticleDef->GetPDGLifeTime();
  G4double aMass = aParticle->GetMass();

#ifdef G4VERBOSE
  if (GetVerboseLevel()>2) {
    G4cout << "G4RadioactiveDecay::GetMeanFreePath() "<< G4endl;
    G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
    G4cout << "Mass:" << aMass/GeV <<"[GeV]";
    G4cout << "c*Tau:" << aCtau/m <<"[m]" <<G4endl;
  }
#endif

  // check if the particle is stable?
  if (aParticleDef->GetPDGStable()) {
    pathlength = DBL_MAX;

  } else if (aCtau < 0.0) {
    pathlength =  DBL_MAX;

    //check if the particle has very short life time ?
  } else if (aCtau < DBL_MIN) {
    pathlength =  DBL_MIN;

    //check if zero mass
  } else if (aMass <  DBL_MIN)  {
    pathlength =  DBL_MAX;
#ifdef G4VERBOSE
    if (GetVerboseLevel()>1) {
      G4cerr << " Zero Mass particle " << G4endl;
    }
#endif
  } else {
    //calculate the mean free path
    // by using normalized kinetic energy (= Ekin/mass)
    G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
    if ( rKineticEnergy > HighestValue) {
      // beta >> 1
      pathlength = ( rKineticEnergy + 1.0)* aCtau;
    } else if ( rKineticEnergy < DBL_MIN ) {
      // too slow particle
#ifdef G4VERBOSE
      if (GetVerboseLevel()>2) {
        G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
        G4cout << aParticleDef->GetParticleName() << G4endl;
        G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV 
               <<"[GeV]";
      }
#endif
      pathlength = DBL_MIN;
    } else {
      // beta << 1
      pathlength = aCtau*(aParticle->GetTotalMomentum())/aMass;
    }
  }
#ifdef G4VERBOSE
  if (GetVerboseLevel()>1) {
    G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
  }
#endif
  return  pathlength;
}

G4double G4RadioactiveDecay::GetMeanLifeTime	(	const G4Track &	theTrack,
		G4ForceCondition *	condition
	)			[protected, virtual]

Implements G4VRestDiscreteProcess.

Definition at line 551 of file G4RadioactiveDecay.cc.

References DBL_MAX, G4cout, G4endl, G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4ParticleDefinition::GetPDGLifeTime(), G4ParticleDefinition::GetPDGStable(), GetVerboseLevel(), and ns.

{
  // For varience reduction implementation the time is set to 0 so as to 
  // force the particle to decay immediately.
  // in analogueMC mode it return the particles meanlife.
  // 
  G4double meanlife = 0.;
  if (AnalogueMC) {
    const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
    G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
    G4double theLife = theParticleDef->GetPDGLifeTime();

#ifdef G4VERBOSE
    if (GetVerboseLevel()>2)
      {
        G4cout <<"G4RadioactiveDecay::GetMeanLifeTime() " <<G4endl;
        G4cout <<"KineticEnergy:" <<theParticle->GetKineticEnergy()/GeV <<"[GeV]";
        G4cout <<"Mass:" <<theParticle->GetMass()/GeV <<"[GeV]"; 
        G4cout <<"Life time: " <<theLife/ns <<"[ns]" << G4endl;
      }
#endif
    if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
    else if (theLife < 0.0) {meanlife = DBL_MAX;}
    else {meanlife = theLife;}
    // set the meanlife to zero for excited istopes which is not in the RDM database
    if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
  }
#ifdef G4VERBOSE
  if (GetVerboseLevel()>1)
    {G4cout <<"mean life time: " <<meanlife/s <<"[s]" <<G4endl;}
#endif

  return  meanlife;
}

G4NucleusLimits G4RadioactiveDecay::GetNucleusLimits

(

)

const [inline]

Definition at line 179 of file G4RadioactiveDecay.hh.

      {return theNucleusLimits;}

G4int G4RadioactiveDecay::GetSplitNuclei

(

)

[inline]

Definition at line 210 of file G4RadioactiveDecay.hh.

{return NSplit;}

G4double G4RadioactiveDecay::GetTaoTime	(	const	G4double,
		const	G4double
	)			[protected]

Definition at line 371 of file G4RadioactiveDecay.cc.

References G4cout, G4endl, and GetVerboseLevel().

Referenced by DecayIt().

{
  long double taotime =0.L;
  G4int nbin;
  if ( t > SBin[NSourceBin]) {
    nbin  = NSourceBin;}
  else {
    nbin = 0;
    while (t > SBin[nbin]) nbin++;
    nbin--;}
  long double lt = t ;
  long double ltao = tao;

  if (nbin > 0) { 
    for (G4int i = 0; i < nbin; i++) 
      {
        taotime += (long double)SProfile[i] * (std::exp(-(lt-(long double)SBin[i+1])/ltao)-std::exp(-(lt-(long double)SBin[i])/ltao));
      }
  }
  taotime +=  (long double)SProfile[nbin] * (1.L-std::exp(-(lt-(long double)SBin[nbin])/ltao));
  if (taotime < 0.)  {
    G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
    G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
    G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
    taotime = 0.;
  }
#ifdef G4VERBOSE
  if (GetVerboseLevel()>1)
    {G4cout <<" Tao time: " <<taotime <<G4endl;}
#endif
  return (G4double)taotime ;
}

std::vector<G4RadioactivityTable*> G4RadioactiveDecay::GetTheRadioactivityTables

(

)

[inline]

Definition at line 154 of file G4RadioactiveDecay.hh.

       {return theRadioactivityTables;}

G4int G4RadioactiveDecay::GetVerboseLevel

(

)

const [inline]

Reimplemented from G4VProcess.

Definition at line 171 of file G4RadioactiveDecay.hh.

Referenced by AddDecayRateTable(), ChooseCollimationDirection(), CollimateDecay(), CollimateDecayProduct(), DecayIt(), DeselectAllVolumes(), DeselectAVolume(), DoDecay(), G4RadioactiveDecay(), GetDecayRateTable(), GetDecayTime(), GetMeanFreePath(), GetMeanLifeTime(), GetTaoTime(), LoadDecayTable(), SelectAllVolumes(), SelectAVolume(), SetDecayBias(), and SetSourceTimeProfile().

{return verboseLevel;}

G4bool G4RadioactiveDecay::IsAnalogueMonteCarlo

(

)

[inline]

Definition at line 194 of file G4RadioactiveDecay.hh.

{return AnalogueMC;}

G4bool G4RadioactiveDecay::IsApplicable

(

const G4ParticleDefinition &

)

[virtual]

Reimplemented from G4VProcess.

Definition at line 214 of file G4RadioactiveDecay.cc.

References G4NucleusLimits::GetAMax(), G4NucleusLimits::GetAMin(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetParticleType(), G4ParticleDefinition::GetPDGLifeTime(), G4NucleusLimits::GetZMax(), and G4NucleusLimits::GetZMin().

Referenced by AddDecayRateTable(), and DecayIt().

{
  // All particles, other than G4Ions, are rejected by default
  if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
  if (aParticle.GetParticleName() == "GenericIon") {
    return true;
  } else if (!(aParticle.GetParticleType() == "nucleus")
             || aParticle.GetPDGLifeTime() < 0. ) {
    return false;
  }

  // Determine whether the nuclide falls into the correct A and Z range

  G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
  G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
  if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
    {return false;}
  else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
    {return false;}
  return true;
}

G4bool G4RadioactiveDecay::IsLoaded

(

const G4ParticleDefinition &

)

Definition at line 237 of file G4RadioactiveDecay.cc.

References G4ParticleDefinition::GetParticleName().

Referenced by AddDecayRateTable(), and DecayIt().

{
  // Check whether the radioactive decay data on the ion have already been
  // loaded

  return std::binary_search(LoadedNuclei.begin(),
                            LoadedNuclei.end(),
                            aParticle.GetParticleName());
}

G4bool G4RadioactiveDecay::IsRateTableReady

(

const G4ParticleDefinition &

)

Definition at line 336 of file G4RadioactiveDecay.cc.

References G4ParticleDefinition::GetParticleName().

Referenced by DecayIt().

{
  // Check whether the radioactive decay rates table for the ion has already
  // been calculated.
  G4String aParticleName = aParticle.GetParticleName();
  for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
    if (theDecayRateTableVector[i].GetIonName() == aParticleName)
      return true;
  }
  return false;
}

G4DecayTable * G4RadioactiveDecay::LoadDecayTable

(

G4ParticleDefinition &

theParentNucleus

)

Definition at line 687 of file G4RadioactiveDecay.cc.

References allowed, Alpha, BetaMinus, BetaPlus, G4DecayTable::DumpInfo(), G4DecayTable::entries(), ERROR, FatalException, G4BetaDecayCorrections::FermiFunction(), G4cerr, G4cout, G4endl, G4Exception(), G4VDecayChannel::GetBR(), G4DecayTable::GetDecayChannel(), G4NuclearDecayChannel::GetDecayMode(), G4ParticleDefinition::GetParticleName(), GetVerboseLevel(), G4DecayTable::Insert(), IT, KshellEC, LshellEC, MshellEC, G4NuclearDecayChannel::SetARM(), G4VDecayChannel::SetBR(), G4NuclearDecayChannel::SetHLThreshold(), G4NuclearDecayChannel::SetICM(), G4BetaDecayCorrections::ShapeFactor(), SpFission, and G4String::strip().

Referenced by AddDecayRateTable(), and DecayIt().

{
  // Create and initialise variables used in the method.
  G4DecayTable* theDecayTable = new G4DecayTable();

  // Generate input data file name using Z and A of the parent nucleus
  // file containing radioactive decay data.
  G4int A    = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
  G4int Z    = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
  G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();

  //Check if data have been provided by the user
  G4String file= theUserRadioactiveDataFiles[1000*A+Z];

  if (file =="") {
    if (!getenv("G4RADIOACTIVEDATA") ) {
      G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files."
             << G4endl;
      throw G4HadronicException(__FILE__, __LINE__,
                              "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
    }
    G4String dirName = getenv("G4RADIOACTIVEDATA");

    std::ostringstream os;
    os <<dirName <<"/z" <<Z <<".a" <<A <<'\0';
    file = os.str();
  }

  LoadedNuclei.push_back(theParentNucleus.GetParticleName());
  std::sort( LoadedNuclei.begin(), LoadedNuclei.end() );
  // sort needed to allow binary_search

  std::ifstream DecaySchemeFile(file);

  G4bool found(false);
  if (DecaySchemeFile) { 
    // Initialise variables used for reading in radioactive decay data.
    G4int nMode = 7;
    G4bool modeFirstRecord[7];
    G4double modeTotalBR[7] = {0.0};
    G4double modeSumBR[7];
    for (G4int i = 0; i < nMode; i++) {
      modeFirstRecord[i] = true;
      modeSumBR[i] = 0.0;
    }

    G4bool complete(false);
    char inputChars[100]={' '};
    G4String inputLine;
    G4String recordType("");
    G4RadioactiveDecayMode theDecayMode;
    G4double a(0.0);
    G4double b(0.0);
    G4double c(0.0);
    G4BetaDecayType betaType(allowed);
    G4double e0;

    // Loop through each data file record until you identify the decay
    // data relating to the nuclide of concern.

    while (!complete && !DecaySchemeFile.getline(inputChars, 100).eof()) {
      inputLine = inputChars;
      inputLine = inputLine.strip(1);
      if (inputChars[0] != '#' && inputLine.length() != 0) {
        std::istringstream tmpStream(inputLine);

        if (inputChars[0] == 'P') {
          // Nucleus is a parent type.  Check excitation level to see if it
          // matches that of theParentNucleus
          tmpStream >> recordType >> a >> b;
          if (found) {complete = true;}
          else {found = (std::abs(a*keV - E) < levelTolerance);}

        } else if (found) {
          // The right part of the radioactive decay data file has been found.  Search
          // through it to determine the mode of decay of the subsequent records.
          if (inputChars[0] == 'W') {
#ifdef G4VERBOSE
            if (GetVerboseLevel() > 0) {
              // a comment line identified and print out the message
              //
              G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
              G4cout << "   In data file " << file << G4endl;
              G4cout << "   " << inputLine << G4endl;
            }
#endif
          } else {
            tmpStream >> theDecayMode >> a >> b >> c >> betaType;

            // Allowed transitions are the default. Forbidden transitions are
            // indicated in the last column.
            if (inputLine.length() < 80) betaType = allowed;
            a /= 1000.;
            c /= 1000.;

            switch (theDecayMode) {

              case IT:   // Isomeric transition
              {
                G4ITDecayChannel* anITChannel =
                  new G4ITDecayChannel(GetVerboseLevel(),
                                       (const G4Ions*)& theParentNucleus, b);
                anITChannel->SetICM(applyICM);
                anITChannel->SetARM(applyARM);
                anITChannel->SetHLThreshold(halflifethreshold);
                theDecayTable->Insert(anITChannel);
              }
              break;

              case BetaMinus:
              {
                if (modeFirstRecord[1]) {
                  modeFirstRecord[1] = false;
                  modeTotalBR[1] = b;
                } else {
                  if (c > 0.) {
                    e0 = c*MeV/0.511;
                    G4BetaDecayCorrections corrections(Z+1, A);

                    // array to store function shape
                    G4int npti = 100;                           
                    G4double* pdf = new G4double[npti];

                    G4double e;   // Total electron energy in units of electron mass
                    G4double p;   // Electron momentum in units of electron mass 
                    G4double f;   // Spectral shape function value
                    for (G4int ptn = 0; ptn < npti; ptn++) {
                      // Calculate simple phase space spectrum
                      e = 1. + e0*(ptn+0.5)/100.;
                      p = std::sqrt(e*e - 1.);
                      f = p*e*(e0-e+1)*(e0-e+1);

                      // Apply Fermi factor to get allowed shape
                      f *= corrections.FermiFunction(e);

                      // Apply shape factor for forbidden transitions
                      f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
                      pdf[ptn] = f;
                    }

                    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
                    G4BetaMinusDecayChannel *aBetaMinusChannel = new
                    G4BetaMinusDecayChannel(GetVerboseLevel(), &theParentNucleus,
                                            b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
                    aBetaMinusChannel->SetICM(applyICM);
                    aBetaMinusChannel->SetARM(applyARM);
                    aBetaMinusChannel->SetHLThreshold(halflifethreshold);
                    theDecayTable->Insert(aBetaMinusChannel);
                    modeSumBR[1] += b;
                    delete[] pdf;
                  } // c > 0
                } // if not first record
              }
              break;

              case BetaPlus:
              {
                if (modeFirstRecord[2]) {
                  modeFirstRecord[2] = false;
                  modeTotalBR[2] = b;
                } else {
                  e0 = c*MeV/0.511 - 2.;
                  // Need to test e0 for nuclei which have Q < 2Me in their
                  // data files (e.g. z67.a162)
                  if (e0 > 0.) {
                    G4BetaDecayCorrections corrections(1-Z, A);

                    // array to store function shape
                    G4int npti = 100;                           
                    G4double* pdf = new G4double[npti];

                    G4double e;   // Total positron energy in units of electron mass
                    G4double p;   // Positron momentum in units of electron mass
                    G4double f;   // Spectral shape function value
                    for (G4int ptn = 0; ptn < npti; ptn++) {
                      // Calculate simple phase space spectrum
                      e = 1. + e0*(ptn+0.5)/100.;
                      p = std::sqrt(e*e - 1.);
                      f = p*e*(e0-e+1)*(e0-e+1);

                      // Apply Fermi factor to get allowed shape
                      f *= corrections.FermiFunction(e);

                      // Apply shape factor for forbidden transitions
                      f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
                      pdf[ptn] = f;
                    }
                    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
                    G4BetaPlusDecayChannel *aBetaPlusChannel = new 
                    G4BetaPlusDecayChannel(GetVerboseLevel(), &theParentNucleus,
                                           b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
                    aBetaPlusChannel->SetICM(applyICM);
                    aBetaPlusChannel->SetARM(applyARM);
                    aBetaPlusChannel->SetHLThreshold(halflifethreshold);
                    theDecayTable->Insert(aBetaPlusChannel);
                    modeSumBR[2] += b;
                    delete[] pdf;
                  } // if e0 > 0
                } // if not first record
              }
              break;

              case KshellEC:  // K-shell electron capture

                if (modeFirstRecord[3]) {
                  modeFirstRecord[3] = false;
                  modeTotalBR[3] = b;
                } else {
                  G4KshellECDecayChannel* aKECChannel =
                    new G4KshellECDecayChannel(GetVerboseLevel(),
                                               &theParentNucleus,
                                               b, c*MeV, a*MeV);
                  aKECChannel->SetICM(applyICM);
                  aKECChannel->SetARM(applyARM);
                  aKECChannel->SetHLThreshold(halflifethreshold);
                  theDecayTable->Insert(aKECChannel);
                  modeSumBR[3] += b;
                }
                break;

              case LshellEC:  // L-shell electron capture

                if (modeFirstRecord[4]) {
                  modeFirstRecord[4] = false;
                  modeTotalBR[4] = b;
                } else {
                  G4LshellECDecayChannel *aLECChannel = new
                    G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
                                            b, c*MeV, a*MeV);
                  aLECChannel->SetICM(applyICM);
                  aLECChannel->SetARM(applyARM);
                  aLECChannel->SetHLThreshold(halflifethreshold);
                  theDecayTable->Insert(aLECChannel);
                  modeSumBR[4] += b;
                }
                break;

              case MshellEC:  // M-shell electron capture
                              // In this implementation it is added to L-shell case
                if (modeFirstRecord[5]) {
                  modeFirstRecord[5] = false;
                  modeTotalBR[5] = b;
                } else {
                  G4MshellECDecayChannel* aMECChannel =
                    new G4MshellECDecayChannel(GetVerboseLevel(),
                                               &theParentNucleus,
                                               b, c*MeV, a*MeV);
                  aMECChannel->SetICM(applyICM);
                  aMECChannel->SetARM(applyARM);
                  aMECChannel->SetHLThreshold(halflifethreshold);
                  theDecayTable->Insert(aMECChannel);
                  modeSumBR[5] += b;
                }
                break;

              case Alpha:
                  //G4cout<<"Alpha channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;

                  if (modeFirstRecord[6]) {
                  modeFirstRecord[6] = false;
                  modeTotalBR[6] = b;
                } else {
                  G4AlphaDecayChannel* anAlphaChannel =
                     new G4AlphaDecayChannel(GetVerboseLevel(),
                                             &theParentNucleus,
                                             b, c*MeV, a*MeV);
                  anAlphaChannel->SetICM(applyICM);
                  anAlphaChannel->SetARM(applyARM);
                  anAlphaChannel->SetHLThreshold(halflifethreshold);
                  theDecayTable->Insert(anAlphaChannel);
                  modeSumBR[6] += b;
                }
                break;
              case SpFission:
                  //Still needed to be implemented
                  //G4cout<<"Sp fission channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
                  break;
              case ERROR:

              default:
                G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
                            FatalException, "Selected decay mode does not exist");
            }  // switch
          }  // if char == W
        }  // if char == P 
      }  // if char != #
    }  // While

    // Go through the decay table and make sure that the branching ratios are
    // correctly normalised.

    G4VDecayChannel* theChannel = 0;
    G4NuclearDecayChannel* theNuclearDecayChannel = 0;
    G4String mode = "";

    G4double theBR = 0.0;
    for (G4int i = 0; i < theDecayTable->entries(); i++) {
      theChannel = theDecayTable->GetDecayChannel(i);
      theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
      theDecayMode = theNuclearDecayChannel->GetDecayMode();

      if (theDecayMode != IT) {
        theBR = theChannel->GetBR();
        theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
      }
    } 
  }   // if (DecaySchemeFile)   
  DecaySchemeFile.close();

                
  if (!found && E > 0.) {
    // cases where IT cascade for exited isotopes without entry in RDM database
    // Decay mode is isomeric transition.
    //
    G4ITDecayChannel *anITChannel = new G4ITDecayChannel
      (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
    anITChannel->SetICM(applyICM);
    anITChannel->SetARM(applyARM);
    anITChannel->SetHLThreshold(halflifethreshold);
    theDecayTable->Insert(anITChannel);
  } 
  if (!theDecayTable) {
    //
    // There is no radioactive decay data for this nucleus.  Return a null
    // decay table.
    //
    G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
    theDecayTable = 0;
    return theDecayTable;
  }     
  if (theDecayTable && GetVerboseLevel()>1)
    {
      G4cout <<"G4RadioactiveDecay::LoadDecayTable()" << G4endl;
      G4cout << "  No. of  entries: "<< theDecayTable->entries() <<G4endl;
      theDecayTable ->DumpInfo();
    }

  return theDecayTable;
}

void G4RadioactiveDecay::SelectAllVolumes

(

)

Definition at line 300 of file G4RadioactiveDecay.cc.

References G4cout, G4endl, G4LogicalVolumeStore::GetInstance(), G4LogicalVolume::GetName(), and GetVerboseLevel().

Referenced by G4RadioactiveDecay().

{
  G4LogicalVolumeStore *theLogicalVolumes;
  G4LogicalVolume *volume;
  theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
  ValidVolumes.clear();
#ifdef G4VERBOSE
  if (GetVerboseLevel()>0)
    G4cout << " RDM Applies to all Volumes"  << G4endl;
#endif
  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
    volume=(*theLogicalVolumes)[i];
    ValidVolumes.push_back(volume->GetName());    
#ifdef G4VERBOSE
    if (GetVerboseLevel()>0)
      G4cout << "         RDM Applies to Volume "  << volume->GetName() << G4endl;
#endif
  }
  std::sort(ValidVolumes.begin(), ValidVolumes.end());
  // sort needed in order to allow binary_search
  isAllVolumesMode=true;
}

void G4RadioactiveDecay::SelectAVolume

(

const G4String

aVolume

)

Definition at line 248 of file G4RadioactiveDecay.cc.

References G4cerr, G4cout, G4endl, G4LogicalVolumeStore::GetInstance(), G4LogicalVolume::GetName(), and GetVerboseLevel().

{
  G4LogicalVolumeStore* theLogicalVolumes;
  G4LogicalVolume *volume;
  theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
  for (size_t i = 0; i < theLogicalVolumes->size(); i++){
    volume=(*theLogicalVolumes)[i];
    if (volume->GetName() == aVolume) {
      ValidVolumes.push_back(aVolume);
      std::sort(ValidVolumes.begin(), ValidVolumes.end());
      // sort need for performing binary_search
#ifdef G4VERBOSE
      if (GetVerboseLevel()>0)
        G4cout << " RDM Applies to : " << aVolume << G4endl; 
#endif
    }else if(i ==  theLogicalVolumes->size())
      {
        G4cerr << "SelectAVolume: "<<aVolume << " is not a valid logical volume name"<< G4endl; 
      }
  }
}

void G4RadioactiveDecay::SetAnalogueMonteCarlo

(

G4bool

r

)

[inline]

Definition at line 184 of file G4RadioactiveDecay.hh.

Referenced by SetBRBias(), SetDecayBias(), SetSourceTimeProfile(), and SetSplitNuclei().

                                                  { 
      AnalogueMC  = r; 
      if (!AnalogueMC) halflifethreshold = 1e-6*CLHEP::s;
    }

void G4RadioactiveDecay::SetARM

(

G4bool

arm

)

[inline]

Definition at line 126 of file G4RadioactiveDecay.hh.

{applyARM = arm;}

void G4RadioactiveDecay::SetBRBias

(

G4bool

r

)

[inline]

Definition at line 198 of file G4RadioactiveDecay.hh.

References SetAnalogueMonteCarlo().

                                     {
      BRBias = r;
      SetAnalogueMonteCarlo(0);
     }

void G4RadioactiveDecay::SetDecayBias

(

G4String

filename

)

Definition at line 1407 of file G4RadioactiveDecay.cc.

References FatalException, G4cout, G4endl, G4Exception(), GetVerboseLevel(), and SetAnalogueMonteCarlo().

{
  
  std::ifstream infile(filename, std::ios::in);
  if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_003",
                           FatalException, "Unable to open bias data file" );

  G4double bin, flux;
  G4int dWindows = 0;
  G4int i ;

  theRadioactivityTables.clear();

  NDecayBin = -1;
  while (infile >> bin >> flux ) {
    NDecayBin++;
    if (NDecayBin > 99) {
      G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_004",
                  FatalException, "Input bias file too big (>100 rows)" );
    } else {
      DBin[NDecayBin] = bin * s;
      DProfile[NDecayBin] = flux;
      if (flux > 0.) {
        decayWindows[NDecayBin] = dWindows;
        dWindows++;
        G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
        theRadioactivityTables.push_back(rTable);
      }
    }
  }
  for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
  for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
  // converted to accumulated probabilities

  SetAnalogueMonteCarlo(0);
  infile.close();

#ifdef G4VERBOSE
  if (GetVerboseLevel()>1)
    {G4cout <<" Decay Bias Profile  Nbin = " << NDecayBin <<G4endl;}
#endif
}

void G4RadioactiveDecay::SetDecayCollimation	(	const G4ThreeVector &	theDir,
		G4double	halfAngle = 0.*CLHEP::deg
	)			[inline]

Definition at line 227 of file G4RadioactiveDecay.hh.

References SetDecayDirection(), and SetDecayHalfAngle().

                                                                      {
      SetDecayDirection(theDir);
      SetDecayHalfAngle(halfAngle);
    }

void G4RadioactiveDecay::SetDecayDirection

(

const G4ThreeVector &

theDir

)

[inline]

Definition at line 213 of file G4RadioactiveDecay.hh.

Referenced by SetDecayCollimation().

                                                               {
      forceDecayDirection = theDir.unit();
    }

void G4RadioactiveDecay::SetDecayHalfAngle

(

G4double

halfAngle = 0.*CLHEP::deg

)

[inline]

Definition at line 221 of file G4RadioactiveDecay.hh.

Referenced by SetDecayCollimation().

                                                                  {
      forceDecayHalfAngle = std::min(std::max(0.*CLHEP::deg,halfAngle),180.*CLHEP::deg);
    }

void G4RadioactiveDecay::SetDecayRate	(	G4int	,
		G4int	,
		G4double	,
		G4int	,
		std::vector< G4double >	,
		std::vector< G4double >
	)

Definition at line 1048 of file G4RadioactiveDecay.cc.

References G4RadioactiveDecayRate::SetA(), G4RadioactiveDecayRate::SetDecayRateC(), G4RadioactiveDecayRate::SetE(), G4RadioactiveDecayRate::SetGeneration(), G4RadioactiveDecayRate::SetTaos(), and G4RadioactiveDecayRate::SetZ().

Referenced by AddDecayRateTable().

{ 
  //fill the decay rate vector 
  theDecayRate.SetZ(theZ);
  theDecayRate.SetA(theA);
  theDecayRate.SetE(theE);
  theDecayRate.SetGeneration(theG);
  theDecayRate.SetDecayRateC(theRates);
  theDecayRate.SetTaos(theTaos);
}

void G4RadioactiveDecay::SetFBeta

(

G4bool

r

)

[inline]

Definition at line 191 of file G4RadioactiveDecay.hh.

{ FBeta  = r; }

void G4RadioactiveDecay::SetHLThreshold

(

G4double

hl

)

[inline]

Definition at line 120 of file G4RadioactiveDecay.hh.

{halflifethreshold = hl;}

void G4RadioactiveDecay::SetICM

(

G4bool

icm

)

[inline]

Definition at line 123 of file G4RadioactiveDecay.hh.

{applyICM = icm;}

void G4RadioactiveDecay::SetNucleusLimits

(

G4NucleusLimits

theNucleusLimits1

)

[inline]

Definition at line 174 of file G4RadioactiveDecay.hh.

      {theNucleusLimits = theNucleusLimits1 ;}

void G4RadioactiveDecay::SetSourceTimeProfile

(

G4String

filename

)

Definition at line 1368 of file G4RadioactiveDecay.cc.

References FatalException, G4cout, G4endl, G4Exception(), GetVerboseLevel(), and SetAnalogueMonteCarlo().

{
  std::ifstream infile ( filename, std::ios::in );
  if (!infile) {
    G4ExceptionDescription ed;
    ed << " Could not open file " << filename << G4endl; 
    G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
                FatalException, ed);
  }

  G4double bin, flux;
  NSourceBin = -1;
  while (infile >> bin >> flux ) {
    NSourceBin++;
    if (NSourceBin > 99) {
      G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
                  FatalException, "Input source time file too big (>100 rows)");

    } else {
      SBin[NSourceBin] = bin * s;
      SProfile[NSourceBin] = flux;
    }
  }
  SetAnalogueMonteCarlo(0);
  infile.close();

#ifdef G4VERBOSE
  if (GetVerboseLevel()>1)
    {G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;}
#endif
}

void G4RadioactiveDecay::SetSplitNuclei

(

G4int

r

)

[inline]

Definition at line 204 of file G4RadioactiveDecay.hh.

References SetAnalogueMonteCarlo().

                                         {
      NSplit = r;
      SetAnalogueMonteCarlo(0);
    }

void G4RadioactiveDecay::SetVerboseLevel

(

G4int

value

)

[inline]

Reimplemented from G4VProcess.

Definition at line 168 of file G4RadioactiveDecay.hh.

{verboseLevel = value;}

---

The documentation for this class was generated from the following files:

• G4RadioactiveDecay.hh
• G4RadioactiveDecay.cc

---