G4RadioactiveDecay Class Reference

#include <G4RadioactiveDecay.hh>

Inheritance diagram for G4RadioactiveDecay:

G4VRestDiscreteProcess G4VProcess

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 ()
G4DecayTableLoadDecayTable (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 G4ThreeVectorGetDecayDirection () 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

G4VParticleChangeDecayIt (const G4Track &theTrack, const G4Step &theStep)
G4DecayProductsDoDecay (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().

00145  : G4VRestDiscreteProcess(processName, fDecay), HighestValue(20.0),
00146    isInitialised(false), forceDecayDirection(0.,0.,0.), 
00147    forceDecayHalfAngle(0.*deg), verboseLevel(0)
00148 {
00149 #ifdef G4VERBOSE
00150   if (GetVerboseLevel()>1) {
00151     G4cout <<"G4RadioactiveDecay constructor    Name: ";
00152     G4cout <<processName << G4endl;   }
00153 #endif
00154 
00155   SetProcessSubType(fRadioactiveDecay);
00156 
00157   theRadioactiveDecaymessenger = new G4RadioactiveDecaymessenger(this);
00158   theIsotopeTable              = new G4RIsotopeTable();
00159   pParticleChange              = &fParticleChangeForRadDecay;
00160 
00161   // Now register the Isotope table with G4IonTable.
00162   G4IonTable *theIonTable =
00163     (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
00164   G4VIsotopeTable *aVirtualTable = theIsotopeTable;
00165   theIonTable->RegisterIsotopeTable(aVirtualTable);
00166 
00167   // Reset the contents of the list of nuclei for which decay scheme data
00168   // have been loaded.
00169   LoadedNuclei.clear();
00170 
00171   //
00172   //Reset the list of user define data file
00173   //
00174   theUserRadioactiveDataFiles.clear();
00175 
00176   //
00177   //
00178   // Apply default values.
00179   //
00180   NSourceBin  = 1;
00181   SBin[0]     = 0.* s;
00182   SBin[1]     = 1.* s;
00183   SProfile[0] = 1.;
00184   SProfile[1] = 0.;
00185   NDecayBin   = 1;
00186   DBin[0]     = 0. * s ;
00187   DBin[1]     = 1. * s;
00188   DProfile[0] = 1.;
00189   DProfile[1] = 0.;
00190   decayWindows[0] = 0;
00191   G4RadioactivityTable* rTable = new G4RadioactivityTable() ;
00192   theRadioactivityTables.push_back(rTable);
00193   NSplit      = 1;
00194   AnalogueMC  = true ;
00195   FBeta       = false ;
00196   BRBias      = true ;
00197   applyICM    = true ;
00198   applyARM    = true ;
00199   halflifethreshold = -1.*second;
00200   //
00201   // RDM applies to xall logical volumes as default
00202   isAllVolumesMode=true;
00203   SelectAllVolumes();
00204 }

G4RadioactiveDecay::~G4RadioactiveDecay (  ) 

Definition at line 207 of file G4RadioactiveDecay.cc.

00208 {
00209   delete theRadioactiveDecaymessenger;
00210 }


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

01065 {
01066   // 1) To calculate all the coefficiecies required to derive the
01067   //    radioactivities for all progeny of theParentNucleus
01068   //
01069   // 2) Add the coefficiencies to the decay rate table vector 
01070   //
01071 
01072   //
01073   // Create and initialise variables used in the method.
01074   //
01075   theDecayRateVector.clear();
01076 
01077   G4int nGeneration = 0;
01078   std::vector<G4double> rates;
01079   std::vector<G4double> taos;
01080 
01081   // start rate is -1.
01082   // Eq.4.26 of the Technical Note
01083   rates.push_back(-1.);
01084   //
01085   //
01086   G4int A = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
01087   G4int Z = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
01088   G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
01089   G4double tao = theParentNucleus.GetPDGLifeTime();
01090   if (tao < 0.) tao = 1e-100;
01091   taos.push_back(tao);
01092   G4int nEntry = 0;
01093 
01094   //fill the decay rate with the intial isotope data
01095   SetDecayRate(Z,A,E,nGeneration,rates,taos);
01096 
01097   // store the decay rate in decay rate vector
01098   theDecayRateVector.push_back(theDecayRate);
01099   nEntry++;
01100 
01101   // now start treating the sencondary generations..
01102 
01103   G4bool stable = false;
01104   G4int i;
01105   G4int j;
01106   G4VDecayChannel* theChannel = 0;
01107   G4NuclearDecayChannel* theNuclearDecayChannel = 0;
01108   G4ITDecayChannel* theITChannel = 0;
01109   G4BetaMinusDecayChannel *theBetaMinusChannel = 0;
01110   G4BetaPlusDecayChannel *theBetaPlusChannel = 0;
01111   G4AlphaDecayChannel *theAlphaChannel = 0;
01112   G4RadioactiveDecayMode theDecayMode;
01113   G4double theBR = 0.0;
01114   G4int AP = 0;
01115   G4int ZP = 0;
01116   G4int AD = 0;
01117   G4int ZD = 0;
01118   G4double EP = 0.;
01119   std::vector<G4double> TP;
01120   std::vector<G4double> RP;
01121   G4ParticleDefinition *theDaughterNucleus;
01122   G4double daughterExcitation;
01123   G4ParticleDefinition *aParentNucleus;
01124   G4IonTable* theIonTable;
01125   G4DecayTable *aTempDecayTable;
01126   G4double theRate;
01127   G4double TaoPlus;
01128   G4int nS = 0;
01129   G4int nT = nEntry;
01130   G4double brs[7];
01131   //
01132   theIonTable =
01133     (G4IonTable*)(G4ParticleTable::GetParticleTable()->GetIonTable());
01134 
01135   while (!stable) {
01136     nGeneration++;
01137     for (j = nS; j< nT; j++) {
01138       ZP = theDecayRateVector[j].GetZ();
01139       AP = theDecayRateVector[j].GetA();
01140       EP = theDecayRateVector[j].GetE();
01141       RP = theDecayRateVector[j].GetDecayRateC();
01142       TP = theDecayRateVector[j].GetTaos();      
01143       if (GetVerboseLevel()>0){
01144         G4cout <<"G4RadioactiveDecay::AddDecayRateTable : "
01145                << " daughters of ("<< ZP <<", "<<AP<<", "
01146                << EP <<") "
01147                << " are being calculated. "       
01148                <<" generation = "
01149                << nGeneration << G4endl;
01150       }
01151       aParentNucleus = theIonTable->GetIon(ZP,AP,EP);
01152       if (!IsLoaded(*aParentNucleus)){
01153         aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
01154       }
01155         
01156       G4DecayTable* theDecayTable = new G4DecayTable();
01157       aTempDecayTable = aParentNucleus->GetDecayTable();
01158       for (i=0; i< 7; i++) brs[i] = 0.0;
01159 
01160       //
01161       // Go through the decay table and to combine the same decay channels
01162       //
01163       for (i=0; i<aTempDecayTable->entries(); i++) {
01164         theChannel             = aTempDecayTable->GetDecayChannel(i);
01165         theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
01166         theDecayMode           = theNuclearDecayChannel->GetDecayMode();
01167         daughterExcitation = theNuclearDecayChannel->GetDaughterExcitation ();
01168         theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus () ;
01169         AD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
01170         ZD = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();  
01171         G4NuclearLevelManager* levelManager = 
01172            G4NuclearLevelStore::GetInstance()->GetManager(ZD, AD);
01173         if (levelManager->NumberOfLevels() ) {
01174           const G4NuclearLevel* level =
01175               levelManager->NearestLevel (daughterExcitation);
01176 
01177           if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
01178             // Level half-life is in ns and the threshold is set to 1 micros by default, user can set it via the UI command
01179             if (level->HalfLife()*ns >= halflifethreshold ){    
01180               // save the metastable nucleus 
01181               theDecayTable->Insert(theChannel);
01182             } 
01183             else{
01184               brs[theDecayMode] += theChannel->GetBR();
01185             }
01186           }
01187           else {
01188             brs[theDecayMode] += theChannel->GetBR();
01189           }
01190         }
01191         else{
01192           brs[theDecayMode] += theChannel->GetBR();
01193         }
01194       }     
01195       brs[2] = brs[2]+brs[3]+brs[4]+brs[5];
01196       brs[3] = brs[4] =brs[5] =  0.0;
01197       for (i= 0; i<7; i++){
01198         if (brs[i] > 0.) {
01199           switch ( i ) {
01200           case 0:
01201             //
01202             //
01203             // Decay mode is isomeric transition.
01204             //
01205 
01206             theITChannel =  new G4ITDecayChannel
01207               (0, (const G4Ions*) aParentNucleus, brs[0]);
01208 
01209             theDecayTable->Insert(theITChannel);
01210             break;
01211 
01212           case 1:
01213             //
01214             //
01215             // Decay mode is beta-.
01216             //
01217             theBetaMinusChannel = new G4BetaMinusDecayChannel (0, aParentNucleus,
01218                                                                brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
01219             theDecayTable->Insert(theBetaMinusChannel);
01220 
01221             break;
01222 
01223           case 2:
01224             //
01225             //
01226             // Decay mode is beta+ + EC.
01227             //
01228             theBetaPlusChannel = new G4BetaPlusDecayChannel (GetVerboseLevel(), aParentNucleus,
01229                                                              brs[2], 0.*MeV, 0.*MeV, 1, false, 0);
01230             theDecayTable->Insert(theBetaPlusChannel);
01231             break;                    
01232 
01233           case 6:
01234             //
01235             //
01236             // Decay mode is alpha.
01237             //
01238             theAlphaChannel = new G4AlphaDecayChannel(GetVerboseLevel(),
01239                                                       aParentNucleus,
01240                                                       brs[6], 0.*MeV, 0.*MeV);
01241             theDecayTable->Insert(theAlphaChannel);
01242             break;
01243 
01244           default:
01245             break;
01246           }
01247         }
01248       }
01249       // 
01250       // loop over all branches in theDecayTable
01251       //
01252       for (i = 0; i < theDecayTable->entries(); i++){
01253         theChannel = theDecayTable->GetDecayChannel(i);
01254         theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
01255         theBR = theChannel->GetBR();
01256         theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
01257         //  first check if the decay of the original nucleus is an IT channel, if true create a new groud-level nucleus
01258         if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1) {
01259           A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
01260           Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
01261           theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
01262         }
01263         if (IsApplicable(*theDaughterNucleus) &&
01264             theBR && 
01265             aParentNucleus != theDaughterNucleus) { 
01266           // need to make sure daugher has decaytable
01267           if (!IsLoaded(*theDaughterNucleus)){
01268             theDaughterNucleus->SetDecayTable(LoadDecayTable(*theDaughterNucleus));
01269           }
01270           if (theDaughterNucleus->GetDecayTable()->entries() ) {
01271             //
01272             A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
01273             Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
01274             E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
01275 
01276             TaoPlus = theDaughterNucleus->GetPDGLifeTime();
01277             //          cout << TaoPlus <<G4endl;
01278 
01279             if (TaoPlus <= 0.)  TaoPlus = 1e-100;
01280 
01281 
01282 
01283             // first set the taos, one simply need to add to the parent ones
01284             taos.clear();
01285             taos = TP;
01286             size_t k;
01287             //check that TaoPlus differs from other taos from at least 1.e5 relative difference
01288             //for (k = 0; k < TP.size(); k++){
01289             //  if (std::abs((TaoPlus-TP[k])/TP[k])<1.e-5 ) TaoPlus=1.00001*TP[k];
01290             //}
01291             taos.push_back(TaoPlus);
01292             // now calculate the coefficiencies
01293             //
01294             // they are in two parts, first the less than n ones
01295             // Eq 4.24 of the TN
01296             rates.clear();
01297             long double ta1,ta2;
01298             ta2 = (long double)TaoPlus;
01299             for (k = 0; k < RP.size(); k++){
01300               ta1 = (long double)TP[k];
01301               if (ta1 == ta2) {
01302                 theRate = 1.e100;
01303               }else{
01304                 theRate = ta1/(ta1-ta2);}
01305               theRate = theRate * theBR * RP[k];
01306               rates.push_back(theRate);
01307             }
01308             //
01309             // the sencond part: the n:n coefficiency
01310             // Eq 4.25 of the TN.  Note Yn+1 is zero apart from Y1 which is -1 as treated at line 1013
01311             // 
01312             theRate = 0.;
01313             long double aRate, aRate1;
01314             aRate1 = 0.L;
01315             for (k = 0; k < RP.size(); k++){
01316               ta1 = (long double)TP[k];
01317               if (ta1 == ta2 ) {
01318                 aRate = 1.e100;
01319               }else {
01320                 aRate = ta2/(ta1-ta2);}
01321               aRate = aRate * (long double)(theBR * RP[k]);
01322               aRate1 += aRate;
01323             }
01324             theRate = -aRate1;
01325             rates.push_back(theRate);         
01326             SetDecayRate (Z,A,E,nGeneration,rates,taos);
01327             theDecayRateVector.push_back(theDecayRate);
01328             nEntry++;
01329           }
01330         } // end of testing daughter nucleus
01331       } // end of i loop( the branches) 
01332       //      delete theDecayTable;
01333 
01334     } //end of for j loop
01335     nS = nT;
01336     nT = nEntry;
01337     if (nS == nT) stable = true;
01338   }
01339 
01340   //end of while loop
01341   // the calculation completed here
01342 
01343 
01344   // fill the first part of the decay rate table
01345   // which is the name of the original particle (isotope) 
01346   //
01347   theDecayRateTable.SetIonName(theParentNucleus.GetParticleName()); 
01348   //
01349   //
01350   // now fill the decay table with the newly completed decay rate vector
01351   //
01352 
01353   theDecayRateTable.SetItsRates(theDecayRateVector);
01354 
01355   //
01356   // finally add the decayratetable to the tablevector
01357   //
01358   theDecayRateTableVector.push_back(theDecayRateTable);
01359 }

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

01029 { if (Z<1 || A<2) {
01030         G4cout<<"Z and A not valid!"<<G4endl;
01031   }
01032 
01033   std::ifstream DecaySchemeFile(filename);
01034   if (DecaySchemeFile){
01035         G4int ID_ion=A*1000+Z;
01036         theUserRadioactiveDataFiles[ID_ion]=filename;
01037         theIsotopeTable->AddUserDecayDataFile(Z,A,filename);
01038   }
01039   else {
01040         G4cout<<"The file "<<filename<<" does not exist!"<<G4endl;
01041   }
01042 }

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

00671 {
00672   if(!isInitialised) {
00673     isInitialised = true;
00674     G4VAtomDeexcitation* p = G4LossTableManager::Instance()->AtomDeexcitation();
00675     if(p) { p->InitialiseAtomicDeexcitation(); }
00676   }
00677 }

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

01894                                                                    {
01895   if (origin == forceDecayDirection) return origin;     // Don't do collimation
01896   if (forceDecayHalfAngle == 180.*deg) return origin;
01897 
01898   G4ThreeVector dir = forceDecayDirection;
01899 
01900   // Return direction offset by random throw
01901   if (forceDecayHalfAngle > 0.) {
01902     // Generate uniform direction around central axis
01903     G4double phi = 2.*pi*G4UniformRand();
01904     G4double cosMin = std::cos(forceDecayHalfAngle);
01905     G4double cosTheta = (1.-cosMin)*G4UniformRand() + cosMin;   // [cosMin,1.)
01906     
01907     dir.setPhi(dir.phi()+phi);
01908     dir.setTheta(dir.theta()+std::acos(cosTheta));
01909   }
01910 
01911 #ifdef G4VERBOSE
01912   if (GetVerboseLevel()>1)
01913     G4cout << " ChooseCollimationDirection returns " << dir << G4endl;
01914 #endif
01915 
01916   return dir;
01917 }

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

01852                                                                  {
01853   if (origin == forceDecayDirection) return;    // No collimation requested
01854   if (180.*deg == forceDecayHalfAngle) return;
01855   if (0 == products || 0 == products->entries()) return;
01856 
01857 #ifdef G4VERBOSE
01858   if (GetVerboseLevel()>0) G4cout<<"Begin of CollimateDecay..."<<G4endl;
01859 #endif
01860 
01861   // Particles suitable for directional biasing (for if-blocks below)
01862   static const G4ParticleDefinition* electron = G4Electron::Definition();
01863   static const G4ParticleDefinition* positron = G4Positron::Definition();
01864   static const G4ParticleDefinition* neutron  = G4Neutron::Definition();
01865   static const G4ParticleDefinition* gamma    = G4Gamma::Definition();
01866   static const G4ParticleDefinition* alpha    = G4Alpha::Definition();
01867 
01868   G4ThreeVector newDirection;           // Re-use to avoid memory churn
01869   for (G4int i=0; i<products->entries(); i++) {
01870     G4DynamicParticle* daughter = (*products)[i];
01871     const G4ParticleDefinition* daughterType = daughter->GetParticleDefinition();
01872 
01873     if (daughterType == electron || daughterType == positron ||
01874         daughterType == neutron || daughterType == gamma ||
01875         daughterType == alpha) CollimateDecayProduct(daughter);
01876   }
01877 }

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

01879                                                                           {
01880 #ifdef G4VERBOSE
01881   if (GetVerboseLevel() > 1) {
01882     G4cout << "CollimateDecayProduct for daughter "
01883            << daughter->GetParticleDefinition()->GetParticleName() << G4endl;
01884   }
01885 #endif
01886 
01887   G4ThreeVector collimate = ChooseCollimationDirection();
01888   if (origin != collimate) daughter->SetMomentumDirection(collimate);
01889 }

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

01458 {
01459   // Initialize the G4ParticleChange object. Get the particle details and the
01460   // decay table.
01461 
01462 
01463 
01464   fParticleChangeForRadDecay.Initialize(theTrack);
01465   const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
01466 
01467 
01468   G4ParticleDefinition *theParticleDef = theParticle->GetDefinition();
01469 
01470 
01471   // First check whether RDM applies to the current logical volume
01472   if (!isAllVolumesMode){
01473    if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(),
01474                           theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
01475 #ifdef G4VERBOSE
01476     if (GetVerboseLevel()>0) {
01477       G4cout <<"G4RadioactiveDecay::DecayIt : "
01478              << theTrack.GetVolume()->GetLogicalVolume()->GetName()
01479              << " is not selected for the RDM"<< G4endl;
01480       G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
01481       G4cout << " The Valid volumes are " << G4endl;
01482       for (size_t i = 0; i< ValidVolumes.size(); i++) G4cout << ValidVolumes[i] << G4endl;
01483     }
01484 #endif
01485     fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01486 
01487     // Kill the parent particle.
01488 
01489     fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01490     fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01491     ClearNumberOfInteractionLengthLeft();
01492     return &fParticleChangeForRadDecay;
01493    }
01494   }
01495 
01496   // now check is the particle is valid for RDM
01497 
01498   if (!(IsApplicable(*theParticleDef))) { 
01499     //
01500     // The particle is not a Ion or outside the nucleuslimits for decay
01501     //
01502 #ifdef G4VERBOSE
01503     if (GetVerboseLevel()>0) {
01504       G4cerr <<"G4RadioactiveDecay::DecayIt : "
01505              <<theParticleDef->GetParticleName() 
01506              << " is not a valid nucleus for the RDM"<< G4endl;
01507     }
01508 #endif
01509     fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01510 
01511     //
01512     // Kill the parent particle.
01513     //
01514     fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01515     fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01516     ClearNumberOfInteractionLengthLeft();
01517     return &fParticleChangeForRadDecay;
01518   }
01519 
01520   if (!IsLoaded(*theParticleDef))
01521     theParticleDef->SetDecayTable(LoadDecayTable(*theParticleDef));
01522 
01523   G4DecayTable* theDecayTable = theParticleDef->GetDecayTable();
01524 
01525   if (theDecayTable == 0 || theDecayTable->entries() == 0) {
01526       // There are no data in the decay table.  Set the particle change parameters
01527       // to indicate this.
01528 #ifdef G4VERBOSE
01529       if (GetVerboseLevel()>0)
01530         {
01531           G4cerr <<"G4RadioactiveDecay::DecayIt : decay table not defined  for ";
01532           G4cerr <<theParticleDef->GetParticleName() <<G4endl;
01533         }
01534 #endif
01535       fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01536 
01537       // Kill the parent particle.
01538       fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01539       fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01540       ClearNumberOfInteractionLengthLeft();
01541       return &fParticleChangeForRadDecay;
01542 
01543   } else { 
01544     // Data found.  Try to decay nucleus
01545     G4double energyDeposit = 0.0;
01546     G4double finalGlobalTime = theTrack.GetGlobalTime();
01547     G4double finalLocalTime = theTrack.GetLocalTime();
01548     G4int index;
01549     G4ThreeVector currentPosition;
01550     currentPosition = theTrack.GetPosition();
01551 
01552     // Check whether use Analogue or VR implementation
01553     if (AnalogueMC) {
01554 #ifdef G4VERBOSE
01555       if (GetVerboseLevel() > 0)
01556         G4cout <<"DecayIt:  Analogue MC version " << G4endl;
01557 #endif
01558 
01559       G4DecayProducts* products = DoDecay(*theParticleDef);
01560 
01561       // Check if the product is the same as input and kill the track if
01562       // necessary to prevent infinite loop (11/05/10, F.Lei)
01563       if ( products->entries() == 1) {
01564         fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
01565         fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
01566         fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
01567         ClearNumberOfInteractionLengthLeft();
01568         return &fParticleChangeForRadDecay;
01569       }
01570 
01571       // Get parent particle information and boost the decay products to the
01572       // laboratory frame based on this information.
01573 
01574 
01575       //The Parent Energy used for the boost should be the total energy of
01576       // the nucleus of the parent ion without the energy of the shell electrons
01577       // (correction for bug 1359 by L. Desorgher)
01578       G4double ParentEnergy =  theParticle->GetKineticEnergy()+theParticle->GetParticleDefinition()->GetPDGMass();
01579       G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
01580 
01581 
01582 
01583       if (theTrack.GetTrackStatus() == fStopButAlive) { //this condition seems to be always True, further investigation is needed (L.Desorgher)
01584 
01585         // The particle is decayed at rest.
01586         // since the time is still for rest particle in G4 we need to add the
01587         // additional time lapsed between the particle come to rest and the
01588         // actual decay.  This time is simply sampled with the mean-life of
01589         // the particle.  But we need to protect the case PDGTime < 0.
01590         // (F.Lei 11/05/10)
01591         G4double temptime = -std::log( G4UniformRand())
01592                             *theParticleDef->GetPDGLifeTime();
01593         if (temptime < 0.) temptime = 0.; 
01594         finalGlobalTime += temptime;
01595         finalLocalTime += temptime;
01596         energyDeposit += theParticle->GetKineticEnergy();
01597       }
01598       products->Boost( ParentEnergy, ParentDirection);
01599 
01600 
01601 
01602       // Add products in theParticleChangeForRadDecay.
01603       G4int numberOfSecondaries = products->entries();
01604       fParticleChangeForRadDecay.SetNumberOfSecondaries(numberOfSecondaries);
01605 #ifdef G4VERBOSE
01606       if (GetVerboseLevel()>1) {
01607         G4cout <<"G4RadioactiveDecay::DecayIt : Decay vertex :";
01608         G4cout <<" Time: " <<finalGlobalTime/ns <<"[ns]";
01609         G4cout <<" X:" <<(theTrack.GetPosition()).x() /cm <<"[cm]";
01610         G4cout <<" Y:" <<(theTrack.GetPosition()).y() /cm <<"[cm]";
01611         G4cout <<" Z:" <<(theTrack.GetPosition()).z() /cm <<"[cm]";
01612         G4cout << G4endl;
01613         G4cout <<"G4Decay::DecayIt  : decay products in Lab. Frame" <<G4endl;
01614         products->DumpInfo();
01615         products->IsChecked();
01616       }
01617 #endif
01618       for (index=0; index < numberOfSecondaries; index++) {
01619         G4Track* secondary = new G4Track(products->PopProducts(),
01620                                          finalGlobalTime, currentPosition);
01621         secondary->SetGoodForTrackingFlag();
01622         secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
01623         fParticleChangeForRadDecay.AddSecondary(secondary);
01624       }
01625       delete products;
01626       // end of analogue MC algorithm
01627 
01628     } else {
01629       // Variance Reduction Method
01630 #ifdef G4VERBOSE
01631       if (GetVerboseLevel()>0)
01632         G4cout << "DecayIt: Variance Reduction version " << G4endl;
01633 #endif
01634       if (!IsRateTableReady(*theParticleDef)) {
01635         // if the decayrates are not ready, calculate them and 
01636         // add to the rate table vector 
01637         AddDecayRateTable(*theParticleDef);
01638       }
01639       //retrieve the rates 
01640       GetDecayRateTable(*theParticleDef);
01641 
01642       // declare some of the variables required in the implementation
01643       G4ParticleDefinition* parentNucleus;
01644       G4IonTable* theIonTable;
01645       G4int PZ;
01646       G4int PA;
01647       G4double PE;
01648       std::vector<G4double> PT;
01649       std::vector<G4double> PR;
01650       G4double taotime;
01651       long double decayRate;
01652 
01653       size_t i;
01654       size_t j;
01655       G4int numberOfSecondaries;
01656       G4int totalNumberOfSecondaries = 0;
01657       G4double currentTime = 0.;
01658       G4int ndecaych;
01659       G4DynamicParticle* asecondaryparticle;
01660       std::vector<G4DynamicParticle*> secondaryparticles;
01661       std::vector<G4double> pw;
01662       std::vector<G4double> ptime;
01663       pw.clear();
01664       ptime.clear();
01665 
01666       //now apply the nucleus splitting
01667       for (G4int n = 0; n < NSplit; n++) {
01668         // Get the decay time following the decay probability function 
01669         // suppllied by user  
01670         G4double theDecayTime = GetDecayTime();
01671         G4int nbin = GetDecayTimeBin(theDecayTime);
01672             
01673         // calculate the first part of the weight function  
01674         G4double weight1 = 1.; 
01675         if (nbin == 1) {
01676           weight1 = 1./DProfile[nbin-1] 
01677                     *(DBin[nbin]-DBin[nbin-1])/NSplit;
01678         } else if (nbin > 1) {
01679           weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
01680                     *(DBin[nbin]-DBin[nbin-1])/NSplit;
01681         }
01682 
01683         // it should be calculated in seconds
01684         weight1 /= s ;
01685             
01686         // loop over all the possible secondaries of the nucleus
01687         // the first one is itself.
01688         for (i = 0; i < theDecayRateVector.size(); i++) {
01689           PZ = theDecayRateVector[i].GetZ();
01690           PA = theDecayRateVector[i].GetA();
01691           PE = theDecayRateVector[i].GetE();
01692           PT = theDecayRateVector[i].GetTaos();
01693           PR = theDecayRateVector[i].GetDecayRateC();
01694 
01695           // Calculate the decay rate of the isotope
01696           // decayRate is the radioactivity of isotope (PZ,PA,PE) at the 
01697           // time 'theDecayTime'
01698           // it will be used to calculate the statistical weight of the 
01699           // decay products of this isotope
01700 
01701           //  G4cout <<"PA= "<< PA << " PZ= " << PZ << " PE= "<< PE  <<G4endl;
01702           decayRate = 0.L;
01703           for (j = 0; j < PT.size(); j++) {
01704             taotime = GetTaoTime(theDecayTime,PT[j]);
01705             decayRate -= PR[j] * (long double)taotime;
01706             // Eq.4.23 of of the TN
01707             // note the negative here is required as the rate in the
01708             // equation is defined to be negative, 
01709             // i.e. decay away, but we need positive value here.
01710 
01711             // G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t"
01712             //        << decayRate << G4endl;           
01713           }
01714 
01715           // add the isotope to the radioactivity tables
01716           //  G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl;
01717           //  G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl;
01718           theRadioactivityTables[decayWindows[nbin-1]]->AddIsotope(PZ,PA,PE,weight1*decayRate,theTrack.GetWeight());
01719 
01720           // Now calculate the statistical weight
01721           // One needs to fold the source bias function with the decaytime
01722           // also need to include the track weight! (F.Lei, 28/10/10)
01723           G4double weight = weight1*decayRate*theTrack.GetWeight();
01724 
01725 
01726 
01727           // decay the isotope 
01728           theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
01729           parentNucleus = theIonTable->GetIon(PZ,PA,PE);
01730 
01731           // Create a temprary products buffer.
01732           // Its contents to be transfered to the products at the end of the loop
01733           G4DecayProducts* tempprods = 0;
01734 
01735           // Decide whether to apply branching ratio bias or not             
01736           if (BRBias) {
01737             G4DecayTable* decayTable = parentNucleus->GetDecayTable();
01738             ndecaych = G4int(decayTable->entries()*G4UniformRand());
01739             G4VDecayChannel* theDecayChannel = decayTable->GetDecayChannel(ndecaych);
01740             if (theDecayChannel == 0) {
01741               // Decay channel not found.
01742 #ifdef G4VERBOSE
01743               if (GetVerboseLevel()>0) {
01744                 G4cerr << " G4RadioactiveDecay::DoIt : cannot determine decay channel ";
01745                 G4cerr << " for this nucleus; decay as if no biasing active ";
01746                 G4cerr << G4endl;
01747                 decayTable ->DumpInfo();
01748               }
01749 #endif
01750               tempprods = DoDecay(*parentNucleus);  // DHW 6 Dec 2010 - do decay as if no biasing
01751                                                     //           to avoid deref of temppprods = 0
01752             } else {
01753               // A decay channel has been identified, so execute the DecayIt.
01754               G4double tempmass = parentNucleus->GetPDGMass();
01755               tempprods = theDecayChannel->DecayIt(tempmass);
01756               weight *= (theDecayChannel->GetBR())*(decayTable->entries());
01757             }
01758           } else {
01759             tempprods = DoDecay(*parentNucleus);
01760           }
01761 
01762           // save the secondaries for buffers
01763           numberOfSecondaries = tempprods->entries();
01764           currentTime = finalGlobalTime + theDecayTime;
01765           for (index = 0; index < numberOfSecondaries; index++) {
01766             asecondaryparticle = tempprods->PopProducts();
01767             if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5) {
01768               pw.push_back(weight);
01769               ptime.push_back(currentTime);
01770               secondaryparticles.push_back(asecondaryparticle);
01771             }
01772           }
01773           delete tempprods;
01774 
01775         } // end of i loop
01776       } // end of n loop 
01777 
01778       // now deal with the secondaries in the two stl containers
01779       // and submmit them back to the tracking manager
01780       totalNumberOfSecondaries = pw.size();
01781       fParticleChangeForRadDecay.SetNumberOfSecondaries(totalNumberOfSecondaries);
01782       for (index=0; index < totalNumberOfSecondaries; index++) { 
01783         G4Track* secondary = new G4Track(secondaryparticles[index],
01784                                          ptime[index], currentPosition);
01785         secondary->SetGoodForTrackingFlag();       
01786         secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
01787         secondary->SetWeight(pw[index]);           
01788         fParticleChangeForRadDecay.AddSecondary(secondary); 
01789       }
01790       // make sure the original track is set to stop and its kinematic energy collected
01791       // 
01792       //theTrack.SetTrackStatus(fStopButAlive);
01793       //energyDeposit += theParticle->GetKineticEnergy();
01794 
01795     } // End of Variance Reduction 
01796 
01797     // Kill the parent particle
01798     fParticleChangeForRadDecay.ProposeTrackStatus(fStopAndKill) ;
01799     fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(energyDeposit); 
01800     fParticleChangeForRadDecay.ProposeLocalTime(finalLocalTime);
01801     // Reset NumberOfInteractionLengthLeft.
01802     ClearNumberOfInteractionLengthLeft();
01803 
01804     return &fParticleChangeForRadDecay ;
01805   }
01806 } 

void G4RadioactiveDecay::DeselectAllVolumes (  ) 

Definition at line 324 of file G4RadioactiveDecay.cc.

References G4cout, G4endl, and GetVerboseLevel().

00325 {
00326   ValidVolumes.clear();
00327   isAllVolumesMode=false;
00328 #ifdef G4VERBOSE
00329   if (GetVerboseLevel()>0)
00330     G4cout << " RDM removed from all volumes" << G4endl; 
00331 #endif
00332 }

void G4RadioactiveDecay::DeselectAVolume ( const G4String  aVolume  ) 

Definition at line 271 of file G4RadioactiveDecay.cc.

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

00272 {
00273   G4LogicalVolumeStore *theLogicalVolumes;
00274   G4LogicalVolume *volume;
00275   theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
00276   for (size_t i = 0; i < theLogicalVolumes->size(); i++){
00277     volume=(*theLogicalVolumes)[i];
00278     if (volume->GetName() == aVolume) {
00279       std::vector<G4String>::iterator location;
00280       location = std::find(ValidVolumes.begin(),ValidVolumes.end(),aVolume);
00281       if (location != ValidVolumes.end()) {
00282         ValidVolumes.erase(location);
00283         std::sort(ValidVolumes.begin(), ValidVolumes.end());
00284         isAllVolumesMode =false;
00285       } else {
00286         G4cerr << " DeselectVolume:" << aVolume << " is not in the list"<< G4endl; 
00287       }   
00288 #ifdef G4VERBOSE
00289       if (GetVerboseLevel()>0)
00290         G4cout << " DeselectVolume: " << aVolume << " is removed from list"<<G4endl; 
00291 #endif
00292     } else if (i ==  theLogicalVolumes->size()) {
00293       G4cerr << " DeselectVolume:" << aVolume
00294              << "is not a valid logical volume name"<< G4endl; 
00295     }
00296   }
00297 }

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

01811 {
01812   G4DecayProducts* products = 0;
01813 
01814   // follow the decaytable and generate the secondaries...
01815 #ifdef G4VERBOSE
01816   if (GetVerboseLevel()>0) G4cout<<"Begin of DoDecay..."<<G4endl;
01817 #endif
01818 
01819   G4DecayTable* theDecayTable = theParticleDef.GetDecayTable();
01820 
01821   // Choose a decay channel.
01822 #ifdef G4VERBOSE
01823   if (GetVerboseLevel()>0) G4cout <<"Selecte a channel..."<<G4endl;
01824 #endif
01825 
01826   G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel();
01827   if (theDecayChannel == 0) {
01828     // Decay channel not found.
01829     G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
01830     G4cerr <<G4endl;
01831     theDecayTable ->DumpInfo();
01832   } else {
01833     // A decay channel has been identified, so execute the DecayIt.
01834 #ifdef G4VERBOSE
01835     if (GetVerboseLevel()>1) {
01836       G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
01837       G4cerr <<theDecayChannel <<G4endl;
01838     }
01839 #endif
01840     G4double tempmass = theParticleDef.GetPDGMass();
01841     products = theDecayChannel->DecayIt(tempmass);
01842     // Apply directional bias if requested by user
01843     CollimateDecay(products);
01844   }
01845 
01846   return products;
01847 }

const G4ThreeVector& G4RadioactiveDecay::GetDecayDirection (  )  const [inline]

Definition at line 217 of file G4RadioactiveDecay.hh.

00217                                                           {
00218       return forceDecayDirection; 
00219     }

G4double G4RadioactiveDecay::GetDecayHalfAngle (  )  const [inline]

Definition at line 225 of file G4RadioactiveDecay.hh.

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

00353 {
00354   G4String aParticleName = aParticle.GetParticleName();
00355 
00356   for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
00357     if (theDecayRateTableVector[i].GetIonName() == aParticleName) {
00358       theDecayRateVector = theDecayRateTableVector[i].GetItsRates();
00359     }
00360   }
00361 #ifdef G4VERBOSE
00362   if (GetVerboseLevel() > 0) {
00363     G4cout << "The DecayRate Table for "
00364            << aParticleName << " is selected." <<  G4endl;
00365   }
00366 #endif
00367 }

G4double G4RadioactiveDecay::GetDecayTime (  )  [protected]

Definition at line 522 of file G4RadioactiveDecay.cc.

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

Referenced by DecayIt().

00523 {
00524   G4double decaytime = 0.;
00525   G4double rand = G4UniformRand();
00526   G4int i = 0;
00527   while ( DProfile[i] < rand) i++;
00528   rand = G4UniformRand();
00529   decaytime = DBin[i] + rand*(DBin[i+1]-DBin[i]);
00530 #ifdef G4VERBOSE
00531   if (GetVerboseLevel()>1)
00532     {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;}
00533 #endif
00534   return  decaytime;        
00535 }

G4int G4RadioactiveDecay::GetDecayTimeBin ( const G4double  aDecayTime  )  [protected]

Definition at line 538 of file G4RadioactiveDecay.cc.

Referenced by DecayIt().

00539 {
00540   G4int i = 0;
00541   while ( aDecayTime > DBin[i] ) i++;
00542   return  i;
00543 }

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

00595 {
00596   // get particle
00597   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00598 
00599   // returns the mean free path in GEANT4 internal units
00600   G4double pathlength;
00601   G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00602   G4double aCtau = c_light * aParticleDef->GetPDGLifeTime();
00603   G4double aMass = aParticle->GetMass();
00604 
00605 #ifdef G4VERBOSE
00606   if (GetVerboseLevel()>2) {
00607     G4cout << "G4RadioactiveDecay::GetMeanFreePath() "<< G4endl;
00608     G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
00609     G4cout << "Mass:" << aMass/GeV <<"[GeV]";
00610     G4cout << "c*Tau:" << aCtau/m <<"[m]" <<G4endl;
00611   }
00612 #endif
00613 
00614   // check if the particle is stable?
00615   if (aParticleDef->GetPDGStable()) {
00616     pathlength = DBL_MAX;
00617 
00618   } else if (aCtau < 0.0) {
00619     pathlength =  DBL_MAX;
00620 
00621     //check if the particle has very short life time ?
00622   } else if (aCtau < DBL_MIN) {
00623     pathlength =  DBL_MIN;
00624 
00625     //check if zero mass
00626   } else if (aMass <  DBL_MIN)  {
00627     pathlength =  DBL_MAX;
00628 #ifdef G4VERBOSE
00629     if (GetVerboseLevel()>1) {
00630       G4cerr << " Zero Mass particle " << G4endl;
00631     }
00632 #endif
00633   } else {
00634     //calculate the mean free path
00635     // by using normalized kinetic energy (= Ekin/mass)
00636     G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
00637     if ( rKineticEnergy > HighestValue) {
00638       // beta >> 1
00639       pathlength = ( rKineticEnergy + 1.0)* aCtau;
00640     } else if ( rKineticEnergy < DBL_MIN ) {
00641       // too slow particle
00642 #ifdef G4VERBOSE
00643       if (GetVerboseLevel()>2) {
00644         G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
00645         G4cout << aParticleDef->GetParticleName() << G4endl;
00646         G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV 
00647                <<"[GeV]";
00648       }
00649 #endif
00650       pathlength = DBL_MIN;
00651     } else {
00652       // beta << 1
00653       pathlength = aCtau*(aParticle->GetTotalMomentum())/aMass;
00654     }
00655   }
00656 #ifdef G4VERBOSE
00657   if (GetVerboseLevel()>1) {
00658     G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
00659   }
00660 #endif
00661   return  pathlength;
00662 }

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.

00553 {
00554   // For varience reduction implementation the time is set to 0 so as to 
00555   // force the particle to decay immediately.
00556   // in analogueMC mode it return the particles meanlife.
00557   // 
00558   G4double meanlife = 0.;
00559   if (AnalogueMC) {
00560     const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
00561     G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
00562     G4double theLife = theParticleDef->GetPDGLifeTime();
00563 
00564 #ifdef G4VERBOSE
00565     if (GetVerboseLevel()>2)
00566       {
00567         G4cout <<"G4RadioactiveDecay::GetMeanLifeTime() " <<G4endl;
00568         G4cout <<"KineticEnergy:" <<theParticle->GetKineticEnergy()/GeV <<"[GeV]";
00569         G4cout <<"Mass:" <<theParticle->GetMass()/GeV <<"[GeV]"; 
00570         G4cout <<"Life time: " <<theLife/ns <<"[ns]" << G4endl;
00571       }
00572 #endif
00573     if (theParticleDef->GetPDGStable()) {meanlife = DBL_MAX;}
00574     else if (theLife < 0.0) {meanlife = DBL_MAX;}
00575     else {meanlife = theLife;}
00576     // set the meanlife to zero for excited istopes which is not in the RDM database
00577     if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
00578   }
00579 #ifdef G4VERBOSE
00580   if (GetVerboseLevel()>1)
00581     {G4cout <<"mean life time: " <<meanlife/s <<"[s]" <<G4endl;}
00582 #endif
00583 
00584   return  meanlife;
00585 }

G4NucleusLimits G4RadioactiveDecay::GetNucleusLimits (  )  const [inline]

Definition at line 179 of file G4RadioactiveDecay.hh.

00180       {return theNucleusLimits;}

G4int G4RadioactiveDecay::GetSplitNuclei (  )  [inline]

Definition at line 210 of file G4RadioactiveDecay.hh.

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

00372 {
00373   long double taotime =0.L;
00374   G4int nbin;
00375   if ( t > SBin[NSourceBin]) {
00376     nbin  = NSourceBin;}
00377   else {
00378     nbin = 0;
00379     while (t > SBin[nbin]) nbin++;
00380     nbin--;}
00381   long double lt = t ;
00382   long double ltao = tao;
00383 
00384   if (nbin > 0) { 
00385     for (G4int i = 0; i < nbin; i++) 
00386       {
00387         taotime += (long double)SProfile[i] * (std::exp(-(lt-(long double)SBin[i+1])/ltao)-std::exp(-(lt-(long double)SBin[i])/ltao));
00388       }
00389   }
00390   taotime +=  (long double)SProfile[nbin] * (1.L-std::exp(-(lt-(long double)SBin[nbin])/ltao));
00391   if (taotime < 0.)  {
00392     G4cout <<" Tao time =: " <<taotime << " reset to zero!"<<G4endl;
00393     G4cout <<" t = " << t <<" tao = " <<tao <<G4endl;
00394     G4cout << SBin[nbin] << " " <<SBin[0] << G4endl;
00395     taotime = 0.;
00396   }
00397 #ifdef G4VERBOSE
00398   if (GetVerboseLevel()>1)
00399     {G4cout <<" Tao time: " <<taotime <<G4endl;}
00400 #endif
00401   return (G4double)taotime ;
00402 }

std::vector<G4RadioactivityTable*> G4RadioactiveDecay::GetTheRadioactivityTables (  )  [inline]

Definition at line 154 of file G4RadioactiveDecay.hh.

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

00171 {return verboseLevel;}

G4bool G4RadioactiveDecay::IsAnalogueMonteCarlo (  )  [inline]

Definition at line 194 of file G4RadioactiveDecay.hh.

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

00215 {
00216   // All particles, other than G4Ions, are rejected by default
00217   if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
00218   if (aParticle.GetParticleName() == "GenericIon") {
00219     return true;
00220   } else if (!(aParticle.GetParticleType() == "nucleus")
00221              || aParticle.GetPDGLifeTime() < 0. ) {
00222     return false;
00223   }
00224 
00225   // Determine whether the nuclide falls into the correct A and Z range
00226 
00227   G4int A = ((const G4Ions*) (&aParticle))->GetAtomicMass();
00228   G4int Z = ((const G4Ions*) (&aParticle))->GetAtomicNumber();
00229   if (A > theNucleusLimits.GetAMax() || A < theNucleusLimits.GetAMin())
00230     {return false;}
00231   else if (Z > theNucleusLimits.GetZMax() || Z < theNucleusLimits.GetZMin())
00232     {return false;}
00233   return true;
00234 }

G4bool G4RadioactiveDecay::IsLoaded ( const G4ParticleDefinition  ) 

Definition at line 237 of file G4RadioactiveDecay.cc.

References G4ParticleDefinition::GetParticleName().

Referenced by AddDecayRateTable(), and DecayIt().

00238 {
00239   // Check whether the radioactive decay data on the ion have already been
00240   // loaded
00241 
00242   return std::binary_search(LoadedNuclei.begin(),
00243                             LoadedNuclei.end(),
00244                             aParticle.GetParticleName());
00245 }

G4bool G4RadioactiveDecay::IsRateTableReady ( const G4ParticleDefinition  ) 

Definition at line 336 of file G4RadioactiveDecay.cc.

References G4ParticleDefinition::GetParticleName().

Referenced by DecayIt().

00337 {
00338   // Check whether the radioactive decay rates table for the ion has already
00339   // been calculated.
00340   G4String aParticleName = aParticle.GetParticleName();
00341   for (size_t i = 0; i < theDecayRateTableVector.size(); i++) {
00342     if (theDecayRateTableVector[i].GetIonName() == aParticleName)
00343       return true;
00344   }
00345   return false;
00346 }

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

00688 {
00689   // Create and initialise variables used in the method.
00690   G4DecayTable* theDecayTable = new G4DecayTable();
00691 
00692   // Generate input data file name using Z and A of the parent nucleus
00693   // file containing radioactive decay data.
00694   G4int A    = ((const G4Ions*)(&theParentNucleus))->GetAtomicMass();
00695   G4int Z    = ((const G4Ions*)(&theParentNucleus))->GetAtomicNumber();
00696   G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
00697 
00698   //Check if data have been provided by the user
00699   G4String file= theUserRadioactiveDataFiles[1000*A+Z];
00700 
00701   if (file =="") {
00702     if (!getenv("G4RADIOACTIVEDATA") ) {
00703       G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files."
00704              << G4endl;
00705       throw G4HadronicException(__FILE__, __LINE__,
00706                               "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
00707     }
00708     G4String dirName = getenv("G4RADIOACTIVEDATA");
00709 
00710     std::ostringstream os;
00711     os <<dirName <<"/z" <<Z <<".a" <<A <<'\0';
00712     file = os.str();
00713   }
00714 
00715   LoadedNuclei.push_back(theParentNucleus.GetParticleName());
00716   std::sort( LoadedNuclei.begin(), LoadedNuclei.end() );
00717   // sort needed to allow binary_search
00718 
00719   std::ifstream DecaySchemeFile(file);
00720 
00721   G4bool found(false);
00722   if (DecaySchemeFile) { 
00723     // Initialise variables used for reading in radioactive decay data.
00724     G4int nMode = 7;
00725     G4bool modeFirstRecord[7];
00726     G4double modeTotalBR[7] = {0.0};
00727     G4double modeSumBR[7];
00728     for (G4int i = 0; i < nMode; i++) {
00729       modeFirstRecord[i] = true;
00730       modeSumBR[i] = 0.0;
00731     }
00732 
00733     G4bool complete(false);
00734     char inputChars[100]={' '};
00735     G4String inputLine;
00736     G4String recordType("");
00737     G4RadioactiveDecayMode theDecayMode;
00738     G4double a(0.0);
00739     G4double b(0.0);
00740     G4double c(0.0);
00741     G4BetaDecayType betaType(allowed);
00742     G4double e0;
00743 
00744     // Loop through each data file record until you identify the decay
00745     // data relating to the nuclide of concern.
00746 
00747     while (!complete && !DecaySchemeFile.getline(inputChars, 100).eof()) {
00748       inputLine = inputChars;
00749       inputLine = inputLine.strip(1);
00750       if (inputChars[0] != '#' && inputLine.length() != 0) {
00751         std::istringstream tmpStream(inputLine);
00752 
00753         if (inputChars[0] == 'P') {
00754           // Nucleus is a parent type.  Check excitation level to see if it
00755           // matches that of theParentNucleus
00756           tmpStream >> recordType >> a >> b;
00757           if (found) {complete = true;}
00758           else {found = (std::abs(a*keV - E) < levelTolerance);}
00759 
00760         } else if (found) {
00761           // The right part of the radioactive decay data file has been found.  Search
00762           // through it to determine the mode of decay of the subsequent records.
00763           if (inputChars[0] == 'W') {
00764 #ifdef G4VERBOSE
00765             if (GetVerboseLevel() > 0) {
00766               // a comment line identified and print out the message
00767               //
00768               G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
00769               G4cout << "   In data file " << file << G4endl;
00770               G4cout << "   " << inputLine << G4endl;
00771             }
00772 #endif
00773           } else {
00774             tmpStream >> theDecayMode >> a >> b >> c >> betaType;
00775 
00776             // Allowed transitions are the default. Forbidden transitions are
00777             // indicated in the last column.
00778             if (inputLine.length() < 80) betaType = allowed;
00779             a /= 1000.;
00780             c /= 1000.;
00781 
00782             switch (theDecayMode) {
00783 
00784               case IT:   // Isomeric transition
00785               {
00786                 G4ITDecayChannel* anITChannel =
00787                   new G4ITDecayChannel(GetVerboseLevel(),
00788                                        (const G4Ions*)& theParentNucleus, b);
00789                 anITChannel->SetICM(applyICM);
00790                 anITChannel->SetARM(applyARM);
00791                 anITChannel->SetHLThreshold(halflifethreshold);
00792                 theDecayTable->Insert(anITChannel);
00793               }
00794               break;
00795 
00796               case BetaMinus:
00797               {
00798                 if (modeFirstRecord[1]) {
00799                   modeFirstRecord[1] = false;
00800                   modeTotalBR[1] = b;
00801                 } else {
00802                   if (c > 0.) {
00803                     e0 = c*MeV/0.511;
00804                     G4BetaDecayCorrections corrections(Z+1, A);
00805 
00806                     // array to store function shape
00807                     G4int npti = 100;                           
00808                     G4double* pdf = new G4double[npti];
00809 
00810                     G4double e;   // Total electron energy in units of electron mass
00811                     G4double p;   // Electron momentum in units of electron mass 
00812                     G4double f;   // Spectral shape function value
00813                     for (G4int ptn = 0; ptn < npti; ptn++) {
00814                       // Calculate simple phase space spectrum
00815                       e = 1. + e0*(ptn+0.5)/100.;
00816                       p = std::sqrt(e*e - 1.);
00817                       f = p*e*(e0-e+1)*(e0-e+1);
00818 
00819                       // Apply Fermi factor to get allowed shape
00820                       f *= corrections.FermiFunction(e);
00821 
00822                       // Apply shape factor for forbidden transitions
00823                       f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
00824                       pdf[ptn] = f;
00825                     }
00826 
00827                     RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
00828                     G4BetaMinusDecayChannel *aBetaMinusChannel = new
00829                     G4BetaMinusDecayChannel(GetVerboseLevel(), &theParentNucleus,
00830                                             b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
00831                     aBetaMinusChannel->SetICM(applyICM);
00832                     aBetaMinusChannel->SetARM(applyARM);
00833                     aBetaMinusChannel->SetHLThreshold(halflifethreshold);
00834                     theDecayTable->Insert(aBetaMinusChannel);
00835                     modeSumBR[1] += b;
00836                     delete[] pdf;
00837                   } // c > 0
00838                 } // if not first record
00839               }
00840               break;
00841 
00842               case BetaPlus:
00843               {
00844                 if (modeFirstRecord[2]) {
00845                   modeFirstRecord[2] = false;
00846                   modeTotalBR[2] = b;
00847                 } else {
00848                   e0 = c*MeV/0.511 - 2.;
00849                   // Need to test e0 for nuclei which have Q < 2Me in their
00850                   // data files (e.g. z67.a162)
00851                   if (e0 > 0.) {
00852                     G4BetaDecayCorrections corrections(1-Z, A);
00853 
00854                     // array to store function shape
00855                     G4int npti = 100;                           
00856                     G4double* pdf = new G4double[npti];
00857 
00858                     G4double e;   // Total positron energy in units of electron mass
00859                     G4double p;   // Positron momentum in units of electron mass
00860                     G4double f;   // Spectral shape function value
00861                     for (G4int ptn = 0; ptn < npti; ptn++) {
00862                       // Calculate simple phase space spectrum
00863                       e = 1. + e0*(ptn+0.5)/100.;
00864                       p = std::sqrt(e*e - 1.);
00865                       f = p*e*(e0-e+1)*(e0-e+1);
00866 
00867                       // Apply Fermi factor to get allowed shape
00868                       f *= corrections.FermiFunction(e);
00869 
00870                       // Apply shape factor for forbidden transitions
00871                       f *= corrections.ShapeFactor(betaType, p, e0-e+1.);
00872                       pdf[ptn] = f;
00873                     }
00874                     RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
00875                     G4BetaPlusDecayChannel *aBetaPlusChannel = new 
00876                     G4BetaPlusDecayChannel(GetVerboseLevel(), &theParentNucleus,
00877                                            b, c*MeV, a*MeV, 0, FBeta, aRandomEnergy);
00878                     aBetaPlusChannel->SetICM(applyICM);
00879                     aBetaPlusChannel->SetARM(applyARM);
00880                     aBetaPlusChannel->SetHLThreshold(halflifethreshold);
00881                     theDecayTable->Insert(aBetaPlusChannel);
00882                     modeSumBR[2] += b;
00883                     delete[] pdf;
00884                   } // if e0 > 0
00885                 } // if not first record
00886               }
00887               break;
00888 
00889               case KshellEC:  // K-shell electron capture
00890 
00891                 if (modeFirstRecord[3]) {
00892                   modeFirstRecord[3] = false;
00893                   modeTotalBR[3] = b;
00894                 } else {
00895                   G4KshellECDecayChannel* aKECChannel =
00896                     new G4KshellECDecayChannel(GetVerboseLevel(),
00897                                                &theParentNucleus,
00898                                                b, c*MeV, a*MeV);
00899                   aKECChannel->SetICM(applyICM);
00900                   aKECChannel->SetARM(applyARM);
00901                   aKECChannel->SetHLThreshold(halflifethreshold);
00902                   theDecayTable->Insert(aKECChannel);
00903                   modeSumBR[3] += b;
00904                 }
00905                 break;
00906 
00907               case LshellEC:  // L-shell electron capture
00908 
00909                 if (modeFirstRecord[4]) {
00910                   modeFirstRecord[4] = false;
00911                   modeTotalBR[4] = b;
00912                 } else {
00913                   G4LshellECDecayChannel *aLECChannel = new
00914                     G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
00915                                             b, c*MeV, a*MeV);
00916                   aLECChannel->SetICM(applyICM);
00917                   aLECChannel->SetARM(applyARM);
00918                   aLECChannel->SetHLThreshold(halflifethreshold);
00919                   theDecayTable->Insert(aLECChannel);
00920                   modeSumBR[4] += b;
00921                 }
00922                 break;
00923 
00924               case MshellEC:  // M-shell electron capture
00925                               // In this implementation it is added to L-shell case
00926                 if (modeFirstRecord[5]) {
00927                   modeFirstRecord[5] = false;
00928                   modeTotalBR[5] = b;
00929                 } else {
00930                   G4MshellECDecayChannel* aMECChannel =
00931                     new G4MshellECDecayChannel(GetVerboseLevel(),
00932                                                &theParentNucleus,
00933                                                b, c*MeV, a*MeV);
00934                   aMECChannel->SetICM(applyICM);
00935                   aMECChannel->SetARM(applyARM);
00936                   aMECChannel->SetHLThreshold(halflifethreshold);
00937                   theDecayTable->Insert(aMECChannel);
00938                   modeSumBR[5] += b;
00939                 }
00940                 break;
00941 
00942               case Alpha:
00943                   //G4cout<<"Alpha channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
00944 
00945                   if (modeFirstRecord[6]) {
00946                   modeFirstRecord[6] = false;
00947                   modeTotalBR[6] = b;
00948                 } else {
00949                   G4AlphaDecayChannel* anAlphaChannel =
00950                      new G4AlphaDecayChannel(GetVerboseLevel(),
00951                                              &theParentNucleus,
00952                                              b, c*MeV, a*MeV);
00953                   anAlphaChannel->SetICM(applyICM);
00954                   anAlphaChannel->SetARM(applyARM);
00955                   anAlphaChannel->SetHLThreshold(halflifethreshold);
00956                   theDecayTable->Insert(anAlphaChannel);
00957                   modeSumBR[6] += b;
00958                 }
00959                 break;
00960               case SpFission:
00961                   //Still needed to be implemented
00962                   //G4cout<<"Sp fission channel"<<a<<'\t'<<b<<'\t'<<c<<std::endl;
00963                   break;
00964               case ERROR:
00965 
00966               default:
00967                 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "HAD_RDM_000",
00968                             FatalException, "Selected decay mode does not exist");
00969             }  // switch
00970           }  // if char == W
00971         }  // if char == P 
00972       }  // if char != #
00973     }  // While
00974 
00975     // Go through the decay table and make sure that the branching ratios are
00976     // correctly normalised.
00977 
00978     G4VDecayChannel* theChannel = 0;
00979     G4NuclearDecayChannel* theNuclearDecayChannel = 0;
00980     G4String mode = "";
00981 
00982     G4double theBR = 0.0;
00983     for (G4int i = 0; i < theDecayTable->entries(); i++) {
00984       theChannel = theDecayTable->GetDecayChannel(i);
00985       theNuclearDecayChannel = static_cast<G4NuclearDecayChannel*>(theChannel);
00986       theDecayMode = theNuclearDecayChannel->GetDecayMode();
00987 
00988       if (theDecayMode != IT) {
00989         theBR = theChannel->GetBR();
00990         theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
00991       }
00992     } 
00993   }   // if (DecaySchemeFile)   
00994   DecaySchemeFile.close();
00995 
00996                 
00997   if (!found && E > 0.) {
00998     // cases where IT cascade for exited isotopes without entry in RDM database
00999     // Decay mode is isomeric transition.
01000     //
01001     G4ITDecayChannel *anITChannel = new G4ITDecayChannel
01002       (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
01003     anITChannel->SetICM(applyICM);
01004     anITChannel->SetARM(applyARM);
01005     anITChannel->SetHLThreshold(halflifethreshold);
01006     theDecayTable->Insert(anITChannel);
01007   } 
01008   if (!theDecayTable) {
01009     //
01010     // There is no radioactive decay data for this nucleus.  Return a null
01011     // decay table.
01012     //
01013     G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
01014     theDecayTable = 0;
01015     return theDecayTable;
01016   }     
01017   if (theDecayTable && GetVerboseLevel()>1)
01018     {
01019       G4cout <<"G4RadioactiveDecay::LoadDecayTable()" << G4endl;
01020       G4cout << "  No. of  entries: "<< theDecayTable->entries() <<G4endl;
01021       theDecayTable ->DumpInfo();
01022     }
01023 
01024   return theDecayTable;
01025 }

void G4RadioactiveDecay::SelectAllVolumes (  ) 

Definition at line 300 of file G4RadioactiveDecay.cc.

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

Referenced by G4RadioactiveDecay().

00301 {
00302   G4LogicalVolumeStore *theLogicalVolumes;
00303   G4LogicalVolume *volume;
00304   theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
00305   ValidVolumes.clear();
00306 #ifdef G4VERBOSE
00307   if (GetVerboseLevel()>0)
00308     G4cout << " RDM Applies to all Volumes"  << G4endl;
00309 #endif
00310   for (size_t i = 0; i < theLogicalVolumes->size(); i++){
00311     volume=(*theLogicalVolumes)[i];
00312     ValidVolumes.push_back(volume->GetName());    
00313 #ifdef G4VERBOSE
00314     if (GetVerboseLevel()>0)
00315       G4cout << "         RDM Applies to Volume "  << volume->GetName() << G4endl;
00316 #endif
00317   }
00318   std::sort(ValidVolumes.begin(), ValidVolumes.end());
00319   // sort needed in order to allow binary_search
00320   isAllVolumesMode=true;
00321 }

void G4RadioactiveDecay::SelectAVolume ( const G4String  aVolume  ) 

Definition at line 248 of file G4RadioactiveDecay.cc.

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

00249 {
00250   G4LogicalVolumeStore* theLogicalVolumes;
00251   G4LogicalVolume *volume;
00252   theLogicalVolumes=G4LogicalVolumeStore::GetInstance();
00253   for (size_t i = 0; i < theLogicalVolumes->size(); i++){
00254     volume=(*theLogicalVolumes)[i];
00255     if (volume->GetName() == aVolume) {
00256       ValidVolumes.push_back(aVolume);
00257       std::sort(ValidVolumes.begin(), ValidVolumes.end());
00258       // sort need for performing binary_search
00259 #ifdef G4VERBOSE
00260       if (GetVerboseLevel()>0)
00261         G4cout << " RDM Applies to : " << aVolume << G4endl; 
00262 #endif
00263     }else if(i ==  theLogicalVolumes->size())
00264       {
00265         G4cerr << "SelectAVolume: "<<aVolume << " is not a valid logical volume name"<< G4endl; 
00266       }
00267   }
00268 }

void G4RadioactiveDecay::SetAnalogueMonteCarlo ( G4bool  r  )  [inline]

Definition at line 184 of file G4RadioactiveDecay.hh.

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

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

void G4RadioactiveDecay::SetARM ( G4bool  arm  )  [inline]

Definition at line 126 of file G4RadioactiveDecay.hh.

00126 {applyARM = arm;}

void G4RadioactiveDecay::SetBRBias ( G4bool  r  )  [inline]

Definition at line 198 of file G4RadioactiveDecay.hh.

References SetAnalogueMonteCarlo().

00198                                      {
00199       BRBias = r;
00200       SetAnalogueMonteCarlo(0);
00201      }

void G4RadioactiveDecay::SetDecayBias ( G4String  filename  ) 

Definition at line 1407 of file G4RadioactiveDecay.cc.

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

01408 {
01409   
01410   std::ifstream infile(filename, std::ios::in);
01411   if (!infile) G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_003",
01412                            FatalException, "Unable to open bias data file" );
01413 
01414   G4double bin, flux;
01415   G4int dWindows = 0;
01416   G4int i ;
01417 
01418   theRadioactivityTables.clear();
01419 
01420   NDecayBin = -1;
01421   while (infile >> bin >> flux ) {
01422     NDecayBin++;
01423     if (NDecayBin > 99) {
01424       G4Exception("G4RadioactiveDecay::SetDecayBias()", "HAD_RDM_004",
01425                   FatalException, "Input bias file too big (>100 rows)" );
01426     } else {
01427       DBin[NDecayBin] = bin * s;
01428       DProfile[NDecayBin] = flux;
01429       if (flux > 0.) {
01430         decayWindows[NDecayBin] = dWindows;
01431         dWindows++;
01432         G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
01433         theRadioactivityTables.push_back(rTable);
01434       }
01435     }
01436   }
01437   for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1];
01438   for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin];
01439   // converted to accumulated probabilities
01440 
01441   SetAnalogueMonteCarlo(0);
01442   infile.close();
01443 
01444 #ifdef G4VERBOSE
01445   if (GetVerboseLevel()>1)
01446     {G4cout <<" Decay Bias Profile  Nbin = " << NDecayBin <<G4endl;}
01447 #endif
01448 }

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

Definition at line 227 of file G4RadioactiveDecay.hh.

References SetDecayDirection(), and SetDecayHalfAngle().

00228                                                                       {
00229       SetDecayDirection(theDir);
00230       SetDecayHalfAngle(halfAngle);
00231     }

void G4RadioactiveDecay::SetDecayDirection ( const G4ThreeVector theDir  )  [inline]

Definition at line 213 of file G4RadioactiveDecay.hh.

Referenced by SetDecayCollimation().

00213                                                                {
00214       forceDecayDirection = theDir.unit();
00215     }

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

Definition at line 221 of file G4RadioactiveDecay.hh.

Referenced by SetDecayCollimation().

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

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

01051 { 
01052   //fill the decay rate vector 
01053   theDecayRate.SetZ(theZ);
01054   theDecayRate.SetA(theA);
01055   theDecayRate.SetE(theE);
01056   theDecayRate.SetGeneration(theG);
01057   theDecayRate.SetDecayRateC(theRates);
01058   theDecayRate.SetTaos(theTaos);
01059 }

void G4RadioactiveDecay::SetFBeta ( G4bool  r  )  [inline]

Definition at line 191 of file G4RadioactiveDecay.hh.

00191 { FBeta  = r; }

void G4RadioactiveDecay::SetHLThreshold ( G4double  hl  )  [inline]

Definition at line 120 of file G4RadioactiveDecay.hh.

00120 {halflifethreshold = hl;}

void G4RadioactiveDecay::SetICM ( G4bool  icm  )  [inline]

Definition at line 123 of file G4RadioactiveDecay.hh.

00123 {applyICM = icm;}

void G4RadioactiveDecay::SetNucleusLimits ( G4NucleusLimits  theNucleusLimits1  )  [inline]

Definition at line 174 of file G4RadioactiveDecay.hh.

00175       {theNucleusLimits = theNucleusLimits1 ;}

void G4RadioactiveDecay::SetSourceTimeProfile ( G4String  filename  ) 

Definition at line 1368 of file G4RadioactiveDecay.cc.

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

01369 {
01370   std::ifstream infile ( filename, std::ios::in );
01371   if (!infile) {
01372     G4ExceptionDescription ed;
01373     ed << " Could not open file " << filename << G4endl; 
01374     G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_001",
01375                 FatalException, ed);
01376   }
01377 
01378   G4double bin, flux;
01379   NSourceBin = -1;
01380   while (infile >> bin >> flux ) {
01381     NSourceBin++;
01382     if (NSourceBin > 99) {
01383       G4Exception("G4RadioactiveDecay::SetSourceTimeProfile()", "HAD_RDM_002",
01384                   FatalException, "Input source time file too big (>100 rows)");
01385 
01386     } else {
01387       SBin[NSourceBin] = bin * s;
01388       SProfile[NSourceBin] = flux;
01389     }
01390   }
01391   SetAnalogueMonteCarlo(0);
01392   infile.close();
01393 
01394 #ifdef G4VERBOSE
01395   if (GetVerboseLevel()>1)
01396     {G4cout <<" Source Timeprofile Nbin = " << NSourceBin <<G4endl;}
01397 #endif
01398 }

void G4RadioactiveDecay::SetSplitNuclei ( G4int  r  )  [inline]

Definition at line 204 of file G4RadioactiveDecay.hh.

References SetAnalogueMonteCarlo().

00204                                          {
00205       NSplit = r;
00206       SetAnalogueMonteCarlo(0);
00207     }

void G4RadioactiveDecay::SetVerboseLevel ( G4int  value  )  [inline]

Reimplemented from G4VProcess.

Definition at line 168 of file G4RadioactiveDecay.hh.

00168 {verboseLevel = value;}


The documentation for this class was generated from the following files:
Generated on Mon May 27 17:53:15 2013 for Geant4 by  doxygen 1.4.7