#include <G4RadioactiveDecay.hh>

Inheritance diagram for G4RadioactiveDecay:

Public Member Functions
	G4RadioactiveDecay (const G4String &processName="RadioactiveDecay")

	~G4RadioactiveDecay ()

G4bool	IsApplicable (const G4ParticleDefinition &)

G4DecayTable *	GetDecayTable (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 &)

Public Member Functions inherited from G4VRestDiscreteProcess
	G4VRestDiscreteProcess (const G4String &, G4ProcessType aType=fNotDefined)

	G4VRestDiscreteProcess (G4VRestDiscreteProcess &)

virtual	~G4VRestDiscreteProcess ()

virtual G4double	PostStepGetPhysicalInteractionLength (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)

virtual G4double	AlongStepGetPhysicalInteractionLength (const G4Track &, G4double, G4double, G4double &, G4GPILSelection *)

virtual G4VParticleChange *	AlongStepDoIt (const G4Track &, const G4Step &)

Public Member Functions inherited from G4VProcess
	G4VProcess (const G4String &aName="NoName", G4ProcessType aType=fNotDefined)

	G4VProcess (const G4VProcess &right)

virtual	~G4VProcess ()

G4int	operator== (const G4VProcess &right) const

G4int	operator!= (const G4VProcess &right) const

G4double	GetCurrentInteractionLength () const

void	SetPILfactor (G4double value)

G4double	GetPILfactor () const

G4double	AlongStepGPIL (const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &proposedSafety, G4GPILSelection *selection)

G4double	AtRestGPIL (const G4Track &track, G4ForceCondition *condition)

G4double	PostStepGPIL (const G4Track &track, G4double previousStepSize, G4ForceCondition *condition)

virtual void	PreparePhysicsTable (const G4ParticleDefinition &)

virtual G4bool	StorePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)

virtual G4bool	RetrievePhysicsTable (const G4ParticleDefinition *, const G4String &, G4bool)

const G4String &	GetPhysicsTableFileName (const G4ParticleDefinition *, const G4String &directory, const G4String &tableName, G4bool ascii=false)

const G4String &	GetProcessName () const

G4ProcessType	GetProcessType () const

void	SetProcessType (G4ProcessType)

G4int	GetProcessSubType () const

void	SetProcessSubType (G4int)

virtual void	StartTracking (G4Track *)

virtual void	EndTracking ()

virtual void	SetProcessManager (const G4ProcessManager *)

virtual const G4ProcessManager *	GetProcessManager ()

virtual void	ResetNumberOfInteractionLengthLeft ()

G4double	GetNumberOfInteractionLengthLeft () const

G4double	GetTotalNumberOfInteractionLengthTraversed () const

G4bool	isAtRestDoItIsEnabled () const

G4bool	isAlongStepDoItIsEnabled () const

G4bool	isPostStepDoItIsEnabled () const

virtual void	DumpInfo () const

void	SetVerboseLevel (G4int value)

G4int	GetVerboseLevel () const

virtual void	SetMasterProcess (G4VProcess *masterP)

const G4VProcess *	GetMasterProcess () const

virtual void	BuildWorkerPhysicsTable (const G4ParticleDefinition &part)

virtual void	PrepareWorkerPhysicsTable (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)

Protected Member Functions inherited from G4VProcess
void	SubtractNumberOfInteractionLengthLeft (G4double previousStepSize)

void	ClearNumberOfInteractionLengthLeft ()

Additional Inherited Members
Static Public Member Functions inherited from G4VProcess
static const G4String &	GetProcessTypeName (G4ProcessType)

Protected Attributes inherited from G4VProcess
const G4ProcessManager *	aProcessManager

G4VParticleChange *	pParticleChange

G4ParticleChange	aParticleChange

G4double	theNumberOfInteractionLengthLeft

G4double	currentInteractionLength

G4double	theInitialNumberOfInteractionLength

G4String	theProcessName

G4String	thePhysicsTableFileName

G4ProcessType	theProcessType

G4int	theProcessSubType

G4double	thePILfactor

G4bool	enableAtRestDoIt

G4bool	enableAlongStepDoIt

G4bool	enablePostStepDoIt

G4int	verboseLevel

Detailed Description

Definition at line 82 of file G4RadioactiveDecay.hh.

Constructor & Destructor Documentation

G4RadioactiveDecay::G4RadioactiveDecay ( const G4String & processName = "RadioactiveDecay" )

Definition at line 146 of file G4RadioactiveDecay.cc.

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

  : G4VRestDiscreteProcess(processName, fDecay), isInitialised(false),
    forceDecayDirection(0.,0.,0.), forceDecayHalfAngle(0.*deg), verboseLevel(0)
 {
 #ifdef G4VERBOSE
   if (GetVerboseLevel() > 1) {
     G4cout << "G4RadioactiveDecay constructor: processName = " << processName
            << G4endl;
   }
 #endif
 
   SetProcessSubType(fRadioactiveDecay);
 
   theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
   pParticleChange = &fParticleChangeForRadDecay;
   theIsotopeTable = new G4RIsotopeTable();
 
   // Regsiter the isotope table to the ion table.
   // Although we are touching the ion table, which is shared, we are not
   // adding particles in this operation.  We can therefore do the registration
   // for each instance of the RDM process and do not need to restrict it to the
   // master process.  It's possible that future, more optimized versions of
   // G4IonTable will require to test for master.  
 
   G4IonTable* theIonTable = G4ParticleTable::GetParticleTable()->GetIonTable();
   theIonTable->RegisterIsotopeTable(theIsotopeTable);
 
   // Reset the list of user defined data files
   theUserRadioactiveDataFiles.clear();
 
   // Instantiate the map of decay tables
 #ifdef G4MULTITHREADED
   if(!master_dkmap) master_dkmap = new DecayTableMap;
 #endif
   dkmap = new DecayTableMap;
 
   // 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 = nanosecond;
 
   // RDM applies to xall logical volumes as default
   isAllVolumesMode = true;
   SelectAllVolumes();
 }

G4RadioactiveDecay::~G4RadioactiveDecay ( )

Definition at line 210 of file G4RadioactiveDecay.cc.

 {
   delete theRadioactiveDecaymessenger;
   delete theIsotopeTable;
   for (DecayTableMap::iterator i = dkmap->begin(); i != dkmap->end(); i++) {
     delete i->second;
   }
   dkmap->clear(); 
 }

Member Function Documentation

void G4RadioactiveDecay::AddDecayRateTable ( const G4ParticleDefinition & theParentNucleus )

Definition at line 1073 of file G4RadioactiveDecay.cc.

References G4NuclearLevel::Energy(), G4DecayTable::entries(), G4cout, G4endl, G4VDecayChannel::GetBR(), G4NuclearDecayChannel::GetDaughterExcitation(), G4NuclearDecayChannel::GetDaughterNucleus(), G4DecayTable::GetDecayChannel(), G4NuclearDecayChannel::GetDecayMode(), GetDecayTable(), G4NuclearLevelStore::GetInstance(), G4IonTable::GetIon(), G4ParticleTable::GetIonTable(), G4NuclearLevelStore::GetManager(), G4ParticleDefinition::GetParticleName(), G4ParticleTable::GetParticleTable(), G4ParticleDefinition::GetPDGLifeTime(), GetVerboseLevel(), G4NuclearLevel::HalfLife(), G4DecayTable::Insert(), IsApplicable(), IT, python.hepunit::MeV, G4NuclearLevelManager::NearestLevel(), ns, G4NuclearLevelManager::NumberOfLevels(), SetDecayRate(), 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);
       aTempDecayTable = GetDecayTable(aParentNucleus);
 
       G4DecayTable* theDecayTable = new G4DecayTable();
       for (G4int k = 0; k < 7; k++) brs[k] = 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 daughter has decay table
           aTempDecayTable = GetDecayTable(theDaughterNucleus);
 
           if (aTempDecayTable->entries() ) {
             A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
             Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
             E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
 
             TaoPlus = theDaughterNucleus->GetPDGLifeTime();
             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 1042 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 666 of file G4RadioactiveDecay.cc.

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

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

G4ThreeVector G4RadioactiveDecay::ChooseCollimationDirection ( ) const

protected

Definition at line 1849 of file G4RadioactiveDecay.cc.

References python.hepunit::deg, G4cout, G4endl, G4UniformRand, GetVerboseLevel(), CLHEP::Hep3Vector::phi(), python.hepunit::pi, CLHEP::Hep3Vector::setPhi(), CLHEP::Hep3Vector::setTheta(), and CLHEP::Hep3Vector::theta().

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 1807 of file G4RadioactiveDecay.cc.

References CollimateDecayProduct(), G4Electron::Definition(), G4Positron::Definition(), G4Neutron::Definition(), G4Alpha::Definition(), G4Gamma::Definition(), python.hepunit::deg, G4InuclParticleNames::electron, G4DecayProducts::entries(), G4cout, G4endl, G4DynamicParticle::GetParticleDefinition(), GetVerboseLevel(), neutron, and G4InuclParticleNames::positron.

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 1834 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 1433 of file G4RadioactiveDecay.cc.

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

 {
   // Initialize G4ParticleChange object, get particle details and 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 if particle is valid for RDM
   if (!(IsApplicable(*theParticleDef) ) ) { 
     // Particle is not an ion or is 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;
   }
 
   G4DecayTable* theDecayTable = GetDecayTable(theParticleDef);
 
   if (theDecayTable == 0 || theDecayTable->entries() == 0) {
     // No data in the decay table.  Set 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;
       G4String keyName;
       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 = GetDecayTable(parentNucleus);
 
             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 337 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 282 of file G4RadioactiveDecay.cc.

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

 {
   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 1770 of file G4RadioactiveDecay.cc.

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

Referenced by DecayIt().

 {
   G4DecayProducts* products = 0;
   G4DecayTable* theDecayTable = GetDecayTable(&theParticleDef);
 
   // Choose a decay channel.
 #ifdef G4VERBOSE
   if (GetVerboseLevel() > 0) G4cout << "Select 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;
   } 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.

225 {return forceDecayHalfAngle;}

void G4RadioactiveDecay::GetDecayRateTable ( const G4ParticleDefinition & aParticle )

Definition at line 363 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
 }

G4DecayTable * G4RadioactiveDecay::GetDecayTable ( G4ParticleDefinition * aNucleus )

Definition at line 242 of file G4RadioactiveDecay.cc.

References G4ParticleDefinition::GetParticleName(), and LoadDecayTable().

Referenced by AddDecayRateTable(), DecayIt(), and DoDecay().

 {
   G4String key = aNucleus->GetParticleName();
   DecayTableMap::iterator table_ptr = dkmap->find(key);
 
   G4DecayTable* theDecayTable = 0;
   if (table_ptr == dkmap->end() ) {              // If table not there,     
     theDecayTable = LoadDecayTable(*aNucleus);   // load from file and
     (*dkmap)[key] = theDecayTable;               // store in library 
   } else {
     theDecayTable = table_ptr->second;
   }
 
   return theDecayTable;
 }

G4double G4RadioactiveDecay::GetDecayTime ( )

protected

Definition at line 533 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 549 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
	)

protectedvirtual

Implements G4VRestDiscreteProcess.

Definition at line 605 of file G4RadioactiveDecay.cc.

References python.hepunit::c_light, DBL_MAX, DBL_MIN, G4cout, G4endl, G4Exception(), G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4ParticleDefinition::GetParticleName(), G4ParticleDefinition::GetPDGLifeTime(), G4DynamicParticle::GetTotalMomentum(), GetVerboseLevel(), python.hepunit::GeV, JustWarning, python.hepunit::keV, and python.hepunit::m.

 {
   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
   G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
   G4double tau = aParticleDef->GetPDGLifeTime();
   G4double aMass = aParticle->GetMass();
 
 #ifdef G4VERBOSE
   if (GetVerboseLevel() > 2) {
     G4cout << "G4RadioactiveDecay::GetMeanFreePath() " << G4endl;
     G4cout << "  KineticEnergy: " << aParticle->GetKineticEnergy()/GeV
            << " GeV, Mass: " << aMass/GeV << " GeV, tau: " << tau << " ns "
            << G4endl;
   }
 #endif
   G4double pathlength = DBL_MAX;
   if (tau != -1) {
     // Ion can decay
 
     if (tau < 0.0) {
       G4ExceptionDescription ed;
       ed << "Ion has negative lifetime " << tau
          << " but is not stable.  Setting mean free path to DBL_MAX" << G4endl; 
       G4Exception("G4RadioactiveDecay::GetMeanFreePath()", "HAD_RDM_011",
                    JustWarning, ed);
       pathlength = DBL_MAX;
 
     } else {
       // Calculate mean free path
       G4double betaGamma = aParticle->GetTotalMomentum()/aMass;
       pathlength = c_light*tau*betaGamma;
 
       if (pathlength < DBL_MIN) {
         pathlength = DBL_MIN;
 #ifdef G4VERBOSE
         if (GetVerboseLevel() > 2) {
           G4cout << "G4Decay::GetMeanFreePath: "
                  << aParticleDef->GetParticleName()
                  << " stops, kinetic energy = "
                  << aParticle->GetKineticEnergy()/keV <<" keV " << G4endl;
         }
 #endif
       }
     }
   }
 
 #ifdef G4VERBOSE
   if (GetVerboseLevel() > 1) {
     G4cout << "mean free path: "<< pathlength/m << " m" << G4endl;
   }
 #endif
   return  pathlength;
 }

G4double G4RadioactiveDecay::GetMeanLifeTime	(	const G4Track &	theTrack,
		G4ForceCondition *	condition
	)

protectedvirtual

Implements G4VRestDiscreteProcess.

Definition at line 562 of file G4RadioactiveDecay.cc.

References DBL_MAX, G4cout, G4endl, G4DynamicParticle::GetDefinition(), G4Track::GetDynamicParticle(), G4DynamicParticle::GetKineticEnergy(), G4DynamicParticle::GetMass(), G4ParticleDefinition::GetPDGLifeTime(), G4ParticleDefinition::GetPDGStable(), GetVerboseLevel(), python.hepunit::GeV, 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 particle's mean-life.
 
   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, Mass: " << theParticle->GetMass()/GeV
               << " GeV, 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 meanlife to zero for excited istopes which are 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.

180 {return theNucleusLimits;}

G4int G4RadioactiveDecay::GetSplitNuclei ( )

inline

Definition at line 210 of file G4RadioactiveDecay.hh.

210 {return NSplit;}

G4double G4RadioactiveDecay::GetTaoTime	(	const G4double	t,
		const G4double	tao
	)

protected

Definition at line 382 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 156 of file G4RadioactiveDecay.hh.

Referenced by exrdmAnalysisManager::BeginOfRun(), and exrdmAnalysisManager::EndOfRun().

157 {return theRadioactivityTables;}

G4int G4RadioactiveDecay::GetVerboseLevel ( ) const

inline

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().

171 {return verboseLevel;}

G4bool G4RadioactiveDecay::IsAnalogueMonteCarlo ( )

inline

Definition at line 194 of file G4RadioactiveDecay.hh.

Referenced by exrdmAnalysisManager::BeginOfRun(), and exrdmAnalysisManager::EndOfRun().

194 {return AnalogueMC;}

G4bool G4RadioactiveDecay::IsApplicable ( const G4ParticleDefinition & aParticle )

virtual

Reimplemented from G4VProcess.

Definition at line 221 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::IsRateTableReady ( const G4ParticleDefinition & aParticle )

Definition at line 348 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 689 of file G4RadioactiveDecay.cc.

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

Referenced by GetDecayTable().

 {
   // 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();
   G4int lvl = ((const G4Ions*)(&theParentNucleus))->GetIsomerLevel();
   G4DecayTable* theDecayTable = 0;
 
 #ifdef G4MULTITHREADED
   G4AutoLock lk(&G4RadioactiveDecay::radioactiveDecayMutex);
 
   G4String key = theParentNucleus.GetParticleName();
   DecayTableMap::iterator master_table_ptr = master_dkmap->find(key);
 
   if (master_table_ptr != master_dkmap->end() ) {   // If table is there              
     return master_table_ptr->second;
   }
 #endif
 
   // Create and initialise variables used in the method.
   theDecayTable = new G4DecayTable();
 
   //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();
   }
 
   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);
     G4int levelCounter = 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 {
             found = (levelCounter == lvl);
           } 
           levelCounter++;
 
         } 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;
                     }
 
                     G4RandGeneral* aRandomEnergy = new G4RandGeneral( 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.510999 - 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;
                     }
                     G4RandGeneral* aRandomEnergy = new G4RandGeneral(pdf, npti);  
                     G4BetaPlusDecayChannel* aBetaPlusChannel = new 
                         G4BetaPlusDecayChannel(GetVerboseLevel(),
                                                &theParentNucleus, b,
                                                (c-1.021998)*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:
                   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 RDM_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 && lvl > 0) {
     // Case where IT cascade for excited isotopes has no entries 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;
     theDecayTable->DumpInfo();
   }
 
 #ifdef G4MULTITHREADED
   (*master_dkmap)[key] = theDecayTable;                  // store in master library 
 #endif
   return theDecayTable;
 }

void G4RadioactiveDecay::SelectAllVolumes ( )

Definition at line 313 of file G4RadioactiveDecay.cc.

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

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 259 of file G4RadioactiveDecay.cc.

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

 {
   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 129 of file G4RadioactiveDecay.hh.

Referenced by PhysicsList::AddRadioactiveDecay().

129 {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 1383 of file G4RadioactiveDecay.cc.

References plottest35::bin, 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.

References CLHEP::Hep3Vector::unit().

Referenced by SetDecayCollimation().

                                                                {
       forceDecayDirection = theDir.unit();
     }

void G4RadioactiveDecay::SetDecayHalfAngle ( G4double halfAngle = 0.*CLHEP::deg )

inline

Definition at line 221 of file G4RadioactiveDecay.hh.

References G4INCL::Math::max(), and G4INCL::Math::min().

Referenced by SetDecayCollimation().

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

void G4RadioactiveDecay::SetDecayRate	(	G4int	theZ,
		G4int	theA,
		G4double	theE,
		G4int	theG,
		std::vector< G4double >	theRates,
		std::vector< G4double >	theTaos
	)

Definition at line 1058 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.

191 { FBeta = r; }

void G4RadioactiveDecay::SetHLThreshold ( G4double hl )

inline

Definition at line 123 of file G4RadioactiveDecay.hh.

Referenced by PhysicsList::AddRadioactiveDecay().

123 {halflifethreshold = hl;}

void G4RadioactiveDecay::SetICM ( G4bool icm )

inline

Definition at line 126 of file G4RadioactiveDecay.hh.

Referenced by PhysicsList::AddRadioactiveDecay().

126 {applyICM = icm;}

void G4RadioactiveDecay::SetNucleusLimits ( G4NucleusLimits theNucleusLimits1 )

inline

Definition at line 174 of file G4RadioactiveDecay.hh.

175 {theNucleusLimits = theNucleusLimits1 ;}

void G4RadioactiveDecay::SetSourceTimeProfile ( G4String filename )

Definition at line 1344 of file G4RadioactiveDecay.cc.

References plottest35::bin, 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

Definition at line 168 of file G4RadioactiveDecay.hh.

168 {verboseLevel = value;}

value

const XML_Char int const XML_Char * value

Definition: include/expat.h:331

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

Public Member Functions

Protected Member Functions

Additional Inherited Members

Detailed Description

Constructor & Destructor Documentation

Member Function Documentation