Index: trunk/source/processes/hadronic/models/radioactive_decay/History
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/History	(revision 1340)
+++ trunk/source/processes/hadronic/models/radioactive_decay/History	(revision 1347)
@@ -14,4 +14,38 @@
      * Reverse chronological order (last date on top), please *
      ----------------------------------------------------------
+18 November 2010  F.Lei (radioactive_decay-V09-03-04)
+- src/G4RadioactiveDecay.cc:1593 get rid of the compalition warning with gcc4.5.1
+
+17 November 2010  F.Lei (radioactive_decay-V09-03-03)
+- Completed the implementation of generating the activity table in VR mode
+- Set the default h-l threshold to 1 micros in VR mode.
+- General improvement in VR mode implementation.  
+- G4RadioactiveDecaymessenger.cc: icmCMD,armCMD & hlThreshold are available at all states. 
+
+11 November 2010  F.Lei
+- further updates to G4RadioactivityTable.hh .cc and G4RadioactiveDecay.hh .cc
+
+10 November 2010  Dennis Wright (radioactive_decay-V09-03-02)
+-------------------------------------------------------------
+- G4RadioactiveDecay.cc : replaced incorrect use of "->" with "." for G4Track
+  in DecayIt(const G4Track&, )
+
+- G4RadioactivityTable.cc : add include file for <map> and replace map
+  with std::map in method AddIsotope
+
+29 Oct 2010 F. Lei
+------------------
+- Added G4RadioactivityTable.hh .cc files for tallying the accumulated radioactivitties in VR mode
+- Added in G4RadioactiveDecay.hh:
+	std::vector<G4RadioactivityTable*> GetTheRadioactivityTables() {return theRadioactivityTables;}
+	// this is how the radioactivity tables can be retrieved by the user
+	std::vector<G4RadioactivityTable*>	theRadioactivityTables;
+	std::map<G4int,G4int>		  decayWindows;
+- In G4RadioactiveDecay.cc:
+	starting line 1291: changes need to setup the radioactivity tables
+	line 1624: include the track weight in the weight calculation
+	line 1627: add the rate to the radioactivity tables
+- improved formatting of all the class files
+	
 
 11 Oct 2010 F. Lei (radioactive_decay-V09-03-01)
Index: trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecay.hh
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecay.hh	(revision 1340)
+++ trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecay.hh	(revision 1347)
@@ -61,5 +61,8 @@
 #include "G4RadioactiveDecayRateVector.hh"
 #include "G4RIsotopeTable.hh"
+#include "G4RadioactivityTable.hh"
+
 #include <vector>
+
 
 //class G4UserlimitsForRD;
@@ -74,218 +77,230 @@
 class G4RadioactiveDecay : public G4VRestDiscreteProcess 
 {
-  // class description
-  
-  // Implementation of nuclear radioactive decay process. 
-  // It simulates the decays of radioactive nuclei, which is submitted to RDM in the form
-  // of G4Ions. Half-liffe and decay schemes are hold in the Radioactivity database
-  // All decay products are submitted back to the particle tracking process through 
-  // the G4ParticleChangeForRadDecay object.
-
-  // class description - end 
-  
+	// class description
+
+	// Implementation of nuclear radioactive decay process. 
+	// It simulates the decays of radioactive nuclei, which is submitted to RDM in the form
+	// of G4Ions. Half-liffe and decay schemes are hold in the Radioactivity database
+	// All decay products are submitted back to the particle tracking process through 
+	// the G4ParticleChangeForRadDecay object.
+
+	// class description - end 
+
 public: // with description
 
-  G4RadioactiveDecay (const G4String& processName="RadioactiveDecay");
-  ~G4RadioactiveDecay();
-  // constructor and destructor
-  //
-
-  G4bool IsApplicable(const G4ParticleDefinition &);
-  // Returns true if the specified isotope is
-  //  1) defined as "nucleus" and
-  //  2) it is within theNucleusLimit
-  //
-  G4bool IsLoaded (const G4ParticleDefinition &);
-  // Returns true if the decay table of the specified nuclei is ready.
-  //
-  void SelectAVolume(const G4String aVolume);
-  // Select a logical volume in which RDM applies.
-  //
-  void DeselectAVolume(const G4String aVolume);
-  // remove a logical volume from the RDM applied list
-  //
-  void SelectAllVolumes();
-  // Select all logical volumes for the application of RDM.
-  //
-  void DeselectAllVolumes();
-  // Remove all logical volumes from RDM applications.
-  //
-  void SetDecayBias (G4String filename);
-  //   Sets the decay biasing scheme using the data in "filename"
-  //
-  void SetHLThreshold (G4double hl) {halflifethreshold = hl;}
-  // Set the half-life threshold for isomer production
-  //
-  void SetICM (G4bool icm) {applyICM = icm;}
-  // Enable/disable ICM 
-  //
-  void SetARM (G4bool arm) {applyARM = arm;}
-  // Enable/disable ARM
-  // 
-  void SetSourceTimeProfile (G4String filename) ;
-  //  Sets source exposure function using histograms in "filename"
-  //
-  G4bool IsRateTableReady(const G4ParticleDefinition &);
-  // Returns true if the coefficient and decay time table for all the
-  // descendants of the specified isotope are ready.
-  //
-  // used in VR decay mode only 
-  //
-  void AddDecayRateTable(const G4ParticleDefinition &);
-  // Calculates the coefficient and decay time table for all the descendants
-  // of the specified isotope.  Adds the calculated table to the private data
-  // member "theDecayRateTableVector".
-  //
-  //
-  // used in VR decay mode only 
-  //
-  void GetDecayRateTable(const G4ParticleDefinition &);
-  // Used to retrieve the coefficient and decay time table for all the
-  // descendants of the specified isotope from "theDecayRateTableVector"
-  // and place it in "theDecayRateTable".
-  //
-  //
-  // used in VR decay mode only 
-  //
-  void SetDecayRate(G4int,G4int,G4double, G4int, std::vector<G4double>,
-                    std::vector<G4double>);
-  // Sets "theDecayRate" with data supplied in the arguements.
-  //
-  //  //
-  // used in VR decay mode only 
-  //
-
-  G4DecayTable *LoadDecayTable (G4ParticleDefinition & theParentNucleus);
-  // Load the decay data of isotope theParentNucleus.
-  //
-  inline void  SetVerboseLevel(G4int value) {verboseLevel = value;}
-  // Sets the VerboseLevel which controls duggering display.
-  //
-  inline G4int GetVerboseLevel() const {return verboseLevel;}
-  // Returns the VerboseLevel which controls level of debugging output.
-  //
-  inline void SetNucleusLimits(G4NucleusLimits theNucleusLimits1)
-  {theNucleusLimits = theNucleusLimits1 ;}
-  //  Sets theNucleusLimits which specifies the range of isotopes
-  //  the G4RadioactiveDecay applies.
-  //
-  inline G4NucleusLimits GetNucleusLimits() const
-  {return theNucleusLimits;}
-  //  Returns theNucleusLimits which specifies the range of isotopes
-  //  the G4RadioactiveDecay applies.
-  //
-  inline void SetAnalogueMonteCarlo (G4bool r ) { AnalogueMC  = r; }
-  //   Controls whether G4RadioactiveDecay runs in analogue mode or
-  //   variance reduction mode.
-  inline void SetFBeta (G4bool r ) { FBeta  = r; }
-  //   Controls whether G4RadioactiveDecay uses fast beta simulation mode
-  //
-  inline G4bool IsAnalogueMonteCarlo () {return AnalogueMC;}
-  //   Returns true if the simulation is an analogue Monte Carlo, and false if
-  //   any of the biassing schemes have been selected.
-  //
-  inline void SetBRBias (G4bool r ) { BRBias  = r;
-  SetAnalogueMonteCarlo(0);}
-  //   Sets whether branching ration bias scheme applies.
-  //
-  inline void SetSplitNuclei (G4int r) { NSplit = r; SetAnalogueMonteCarlo(0); }
-  //  Sets the N number for the Nuclei spliting bias scheme
-  //
-  inline G4int GetSplitNuclei () {return NSplit;}
-  //  Returns the N number used for the Nuclei spliting bias scheme
-  //
+	G4RadioactiveDecay (const G4String& processName="RadioactiveDecay");
+	~G4RadioactiveDecay();
+	// constructor and destructor
+	//
+
+	G4bool IsApplicable(const G4ParticleDefinition &);
+	// Returns true if the specified isotope is
+	//  1) defined as "nucleus" and
+	//  2) it is within theNucleusLimit
+	//
+	G4bool IsLoaded (const G4ParticleDefinition &);
+	// Returns true if the decay table of the specified nuclei is ready.
+	//
+	void SelectAVolume(const G4String aVolume);
+	// Select a logical volume in which RDM applies.
+	//
+	void DeselectAVolume(const G4String aVolume);
+	// remove a logical volume from the RDM applied list
+	//
+	void SelectAllVolumes();
+	// Select all logical volumes for the application of RDM.
+	//
+	void DeselectAllVolumes();
+	// Remove all logical volumes from RDM applications.
+	//
+	void SetDecayBias (G4String filename);
+	//   Sets the decay biasing scheme using the data in "filename"
+	//
+	void SetHLThreshold (G4double hl) {halflifethreshold = hl;}
+	// Set the half-life threshold for isomer production
+	//
+	void SetICM (G4bool icm) {applyICM = icm;}
+	// Enable/disable ICM 
+	//
+	void SetARM (G4bool arm) {applyARM = arm;}
+	// Enable/disable ARM
+	// 
+	void SetSourceTimeProfile (G4String filename) ;
+	//  Sets source exposure function using histograms in "filename"
+	//
+	G4bool IsRateTableReady(const G4ParticleDefinition &);
+	// Returns true if the coefficient and decay time table for all the
+	// descendants of the specified isotope are ready.
+	//
+	// used in VR decay mode only 
+	//
+	void AddDecayRateTable(const G4ParticleDefinition &);
+	// Calculates the coefficient and decay time table for all the descendants
+	// of the specified isotope.  Adds the calculated table to the private data
+	// member "theDecayRateTableVector".
+	//
+	//
+	// used in VR decay mode only 
+	//
+	void GetDecayRateTable(const G4ParticleDefinition &);
+	// Used to retrieve the coefficient and decay time table for all the
+	// descendants of the specified isotope from "theDecayRateTableVector"
+	// and place it in "theDecayRateTable".
+	//
+	//
+	// used in VR decay mode only 
+	//
+	void SetDecayRate(G4int,G4int,G4double, G4int, std::vector<G4double>,
+		std::vector<G4double>);
+	// Sets "theDecayRate" with data supplied in the arguements.
+	//
+	//  //
+	// used in VR decay mode only 
+	//
+
+	std::vector<G4RadioactivityTable*> GetTheRadioactivityTables() {return theRadioactivityTables;}
+	// return the vector of G4Radioactivity map
+	// should be used in VR mode only
+
+	G4DecayTable *LoadDecayTable (G4ParticleDefinition & theParentNucleus);
+	// Load the decay data of isotope theParentNucleus.
+	//
+	inline void  SetVerboseLevel(G4int value) {verboseLevel = value;}
+	// Sets the VerboseLevel which controls duggering display.
+	//
+	inline G4int GetVerboseLevel() const {return verboseLevel;}
+	// Returns the VerboseLevel which controls level of debugging output.
+	//
+	inline void SetNucleusLimits(G4NucleusLimits theNucleusLimits1)
+	{theNucleusLimits = theNucleusLimits1 ;}
+	//  Sets theNucleusLimits which specifies the range of isotopes
+	//  the G4RadioactiveDecay applies.
+	//
+	inline G4NucleusLimits GetNucleusLimits() const
+	{return theNucleusLimits;}
+	//  Returns theNucleusLimits which specifies the range of isotopes
+	//  the G4RadioactiveDecay applies.
+	//
+	inline void SetAnalogueMonteCarlo (G4bool r ) { 
+	  AnalogueMC  = r; 
+	  if (!AnalogueMC) halflifethreshold = 1e-6*s;
+	}
+	//   Controls whether G4RadioactiveDecay runs in analogue mode or
+	//   variance reduction mode.
+	inline void SetFBeta (G4bool r ) { FBeta  = r; }
+	//   Controls whether G4RadioactiveDecay uses fast beta simulation mode
+	//
+	inline G4bool IsAnalogueMonteCarlo () {return AnalogueMC;}
+	//   Returns true if the simulation is an analogue Monte Carlo, and false if
+	//   any of the biassing schemes have been selected.
+	//
+	inline void SetBRBias (G4bool r ) { BRBias  = r;
+	SetAnalogueMonteCarlo(0);}
+	//   Sets whether branching ration bias scheme applies.
+	//
+	inline void SetSplitNuclei (G4int r) { NSplit = r; SetAnalogueMonteCarlo(0); }
+	//  Sets the N number for the Nuclei spliting bias scheme
+	//
+	inline G4int GetSplitNuclei () {return NSplit;}
+	//  Returns the N number used for the Nuclei spliting bias scheme
+	//
 public:
-  
-  void BuildPhysicsTable(const G4ParticleDefinition &);
+
+	void BuildPhysicsTable(const G4ParticleDefinition &);
 
 protected:
 
-  G4VParticleChange* DecayIt( const G4Track& theTrack,
-                              const G4Step&  theStep);
-
-  G4DecayProducts* DoDecay(G4ParticleDefinition& theParticleDef);
-
-  G4double GetMeanFreePath(const G4Track& theTrack,
-                           G4double previousStepSize,
-                           G4ForceCondition* condition);
-
-  G4double GetMeanLifeTime(const G4Track& theTrack,
-                           G4ForceCondition* condition);
-
-  G4double GetTaoTime(G4double,G4double);
-
-  G4double GetDecayTime();
-
-  G4int GetDecayTimeBin(const G4double aDecayTime);
+	G4VParticleChange* DecayIt( const G4Track& theTrack,
+		const G4Step&  theStep);
+
+	G4DecayProducts* DoDecay(G4ParticleDefinition& theParticleDef);
+
+	G4double GetMeanFreePath(const G4Track& theTrack,
+		G4double previousStepSize,
+		G4ForceCondition* condition);
+
+	G4double GetMeanLifeTime(const G4Track& theTrack,
+		G4ForceCondition* condition);
+
+	G4double GetTaoTime(G4double,G4double);
+
+	G4double GetDecayTime();
+
+	G4int GetDecayTimeBin(const G4double aDecayTime);
 
 private:
 
-  G4RadioactiveDecay(const G4RadioactiveDecay &right);
-  G4RadioactiveDecay & operator=(const G4RadioactiveDecay &right);
+	G4RadioactiveDecay(const G4RadioactiveDecay &right);
+	G4RadioactiveDecay & operator=(const G4RadioactiveDecay &right);
 
 private:
 
-  G4RadioactiveDecaymessenger  *theRadioactiveDecaymessenger;
-
-  G4PhysicsTable               *aPhysicsTable;
-
-  G4RIsotopeTable              *theIsotopeTable;
-
-  G4NucleusLimits               theNucleusLimits;
-
-  const G4double                HighestBinValue;
-  const G4double                LowestBinValue;
-
-  const G4int                   TotBin;
-
-  G4bool                        AnalogueMC;
-  G4bool                        BRBias;
-  G4bool                        FBeta;
-  G4int                         NSplit;
-
-  G4double                      halflifethreshold;
-  G4bool                        applyICM;
-  G4bool                        applyARM;
-
-  G4int                         NSourceBin;
-  G4double                      SBin[99];
-  G4double                      SProfile[99];
-
-  G4int                         NDecayBin;
-  G4double                      DBin[99];
-  G4double                      DProfile[99];
-
-  std::vector<G4String>         LoadedNuclei;
-  std::vector<G4String>         ValidVolumes;
-
-  G4RadioactiveDecayRate        theDecayRate;
-  G4RadioactiveDecayRates       theDecayRateVector;
-  G4RadioactiveDecayRateVector  theDecayRateTable;
-  G4RadioactiveDecayRateTable   theDecayRateTableVector;
-
-  static const G4double         levelTolerance;
-
-  // Remainder of life time at rest
-  G4double                      fRemainderLifeTime;
-
-  G4int                         verboseLevel;
-
-
-  // ParticleChange for decay process
-  G4ParticleChangeForRadDecay   fParticleChangeForRadDecay;
-
-  inline G4double AtRestGetPhysicalInteractionLength
-    (const G4Track& track, G4ForceCondition* condition)
-  {fRemainderLifeTime = G4VRestDiscreteProcess::
-    AtRestGetPhysicalInteractionLength(track, condition );
-  return fRemainderLifeTime;}
-
-  inline G4VParticleChange* AtRestDoIt
-    (const G4Track& theTrack, const G4Step& theStep)
-  {return DecayIt(theTrack, theStep);}
-
-  inline G4VParticleChange* PostStepDoIt
-    (const G4Track& theTrack, const G4Step& theStep)
-  {return DecayIt(theTrack, theStep);}
+	G4RadioactiveDecaymessenger  *theRadioactiveDecaymessenger;
+
+	G4PhysicsTable               *aPhysicsTable;
+
+	G4RIsotopeTable              *theIsotopeTable;
+
+	G4NucleusLimits               theNucleusLimits;
+
+	const G4double                HighestBinValue;
+	const G4double                LowestBinValue;
+
+	const G4int                   TotBin;
+
+	G4bool                        AnalogueMC;
+	G4bool                        BRBias;
+	G4bool                        FBeta;
+	G4int                         NSplit;
+
+	G4double                      halflifethreshold;
+	G4bool                        applyICM;
+	G4bool                        applyARM;
+
+	G4int                         NSourceBin;
+	G4double                      SBin[99];
+	G4double                      SProfile[99];
+
+	G4int                         NDecayBin;
+	G4double                      DBin[99];
+	G4double                      DProfile[99];
+
+	std::vector<G4String>         LoadedNuclei;
+	std::vector<G4String>         ValidVolumes;
+
+	G4RadioactiveDecayRate        theDecayRate;
+	G4RadioactiveDecayRates       theDecayRateVector;
+        G4RadioactiveDecayRateVector  theDecayRateTable;
+        G4RadioactiveDecayRateTable   theDecayRateTableVector;
+
+	// for the radioactivity tables
+        std::vector<G4RadioactivityTable*>	theRadioactivityTables;
+        G4int		              decayWindows[99];
+
+	//
+	static const G4double		  levelTolerance;
+
+	// Remainder of life time at rest
+	G4double                      fRemainderLifeTime;
+
+	G4int                         verboseLevel;
+
+
+	// ParticleChange for decay process
+	G4ParticleChangeForRadDecay   fParticleChangeForRadDecay;
+
+	inline G4double AtRestGetPhysicalInteractionLength
+		(const G4Track& track, G4ForceCondition* condition)
+	{fRemainderLifeTime = G4VRestDiscreteProcess::
+	AtRestGetPhysicalInteractionLength(track, condition );
+	return fRemainderLifeTime;}
+
+	inline G4VParticleChange* AtRestDoIt
+		(const G4Track& theTrack, const G4Step& theStep)
+	{return DecayIt(theTrack, theStep);}
+
+	inline G4VParticleChange* PostStepDoIt
+		(const G4Track& theTrack, const G4Step& theStep)
+	{return DecayIt(theTrack, theStep);}
 
 };
Index: trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc	(revision 1347)
@@ -121,8 +121,7 @@
 // Constructor
 //
-G4RadioactiveDecay::G4RadioactiveDecay
-  (const G4String& processName)
-  :G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0),
-   LowestBinValue(1.0e-3), TotBin(200), verboseLevel(0)
+G4RadioactiveDecay::G4RadioactiveDecay (const G4String& processName)
+ :G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0),
+  LowestBinValue(1.0e-3), TotBin(200), verboseLevel(0)
 {
 #ifdef G4VERBOSE
@@ -136,5 +135,5 @@
   aPhysicsTable                = 0;
   pParticleChange              = &fParticleChangeForRadDecay;
-  
+
   //
   // Now register the Isotopetable with G4IonTable.
@@ -156,12 +155,15 @@
   NSourceBin  = 1;
   SBin[0]     = 0.* s;
-  SBin[1]     = 1e10 * s;
+  SBin[1]     = 1.* s;
   SProfile[0] = 1.;
-  SProfile[1] = 1.;
+  SProfile[1] = 0.;
   NDecayBin   = 1;
-  DBin[0]     = 9.9e9 * s ;
-  DBin[1]     = 1e10 * s;
+  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 ;
@@ -187,6 +189,6 @@
       delete aPhysicsTable;
     }
-  //  delete theIsotopeTable;
-  delete theRadioactiveDecaymessenger;
+  //    delete theIsotopeTable;
+    delete theRadioactiveDecaymessenger;
 }
 
@@ -197,5 +199,5 @@
 //
 G4bool G4RadioactiveDecay::IsApplicable( const G4ParticleDefinition &
-  aParticle)
+					 aParticle)
 {
   //
@@ -224,5 +226,5 @@
 //
 G4bool G4RadioactiveDecay::IsLoaded(const G4ParticleDefinition &
-  aParticle)
+				    aParticle)
 {
   //
@@ -232,6 +234,6 @@
   //
   return std::binary_search(LoadedNuclei.begin(),
-		       LoadedNuclei.end(),
-		       aParticle.GetParticleName());
+			    LoadedNuclei.end(),
+			    aParticle.GetParticleName());
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -242,5 +244,5 @@
 void G4RadioactiveDecay::SelectAVolume(const G4String aVolume)
 {
-  
+
   G4LogicalVolumeStore *theLogicalVolumes;
   G4LogicalVolume *volume;
@@ -331,5 +333,5 @@
     G4cout << " RDM removed from all volumes" << G4endl; 
 #endif
-  
+
 }
 
@@ -340,5 +342,5 @@
 //
 G4bool G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition &
-  aParticle)
+					    aParticle)
 {
   //
@@ -351,5 +353,5 @@
     {
       if (theDecayRateTableVector[i].GetIonName() == aParticleName)
-        return true ;
+	return true ;
     }
   return false;
@@ -363,5 +365,5 @@
 //
 void G4RadioactiveDecay::GetDecayRateTable(const G4ParticleDefinition &
-  aParticle)
+					   aParticle)
 {
 
@@ -376,9 +378,9 @@
     }
 #ifdef G4VERBOSE
-	if (GetVerboseLevel()>0)
-	  {
-	    G4cout <<"The DecayRate Table for "
-		   << aParticleName << " is selected." <<  G4endl;
-	  }
+  if (GetVerboseLevel()>0)
+    {
+      G4cout <<"The DecayRate Table for "
+	     << aParticleName << " is selected." <<  G4endl;
+    }
 #endif
 }
@@ -388,5 +390,5 @@
 // GetTaoTime
 //
-// to perform the convilution of the sourcetimeprofile function with the 
+// to perform the convolution of the source time profile function with the 
 // decay constants in the decay chain. 
 //
@@ -408,11 +410,16 @@
   }
   taotime +=  SProfile[nbin] * (1-std::exp(-(t-SBin[nbin])/tao));
+  if (taotime < 0.)  {
+    G4cout <<" Tao time: " <<taotime << " reset to zero!"<<G4endl;
+    taotime = 0.;
+  }
+
 #ifdef G4VERBOSE
   if (GetVerboseLevel()>1)
     {G4cout <<" Tao time: " <<taotime <<G4endl;}
 #endif
-  return  taotime ;
-}
- 
+  return taotime ;
+}
+
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -433,5 +440,5 @@
 #ifdef G4VERBOSE
   if (GetVerboseLevel()>1)
-    {G4cout <<" Decay time: " <<decaytime <<"[ns]" <<G4endl;}
+    {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;}
 #endif
   return  decaytime;	    
@@ -447,9 +454,7 @@
 G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime)
 {
-  for (G4int i = 0; i < NDecayBin; i++) 
-    {
-      if ( aDecayTime > DBin[i]) return i+1;	  
-    }
-  return  1;
+  G4int i = 0;
+  while ( aDecayTime > DBin[i] ) i++;
+  return  i;
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -472,5 +477,5 @@
     G4ParticleDefinition* theParticleDef = theParticle->GetDefinition();
     G4double theLife = theParticleDef->GetPDGLifeTime();
-      
+
 #ifdef G4VERBOSE
     if (GetVerboseLevel()>2)
@@ -485,12 +490,12 @@
     else if (theLife < 0.0) {meanlife = DBL_MAX;}
     else {meanlife = theLife;}
-  // set the meanlife to zero for excited istopes which is not in the RDM database
+    // set the meanlife to zero for excited istopes which is not in the RDM database
     if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
   }
 #ifdef G4VERBOSE
   if (GetVerboseLevel()>1)
-    {G4cout <<"mean life time: " <<meanlife <<"[ns]" <<G4endl;}
-#endif
-  
+    {G4cout <<"mean life time: " <<meanlife/s <<"[s]" <<G4endl;}
+#endif
+
   return  meanlife;
 }
@@ -507,8 +512,8 @@
   // constants
   G4bool isOutRange ;
-  
+
   // get particle
   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
-  
+
   // returns the mean free path in GEANT4 internal units
   G4double pathlength;
@@ -516,5 +521,5 @@
   G4double aCtau = c_light * aParticleDef->GetPDGLifeTime();
   G4double aMass = aParticle->GetMass();
-  
+
 #ifdef G4VERBOSE
   if (GetVerboseLevel()>2) {
@@ -525,16 +530,16 @@
   }
 #endif
-  
+
   // check if the particle is stable?
   if (aParticleDef->GetPDGStable()) {
     pathlength = DBL_MAX;
-    
+
   } else if (aCtau < 0.0) {
     pathlength =  DBL_MAX;
-    
+
     //check if the particle has very short life time ?
   } else if (aCtau < DBL_MIN) {
     pathlength =  DBL_MIN;
-    
+
     //check if zero mass
   } else if (aMass <  DBL_MIN)  {
@@ -545,38 +550,38 @@
     }
 #endif
-   } else {
-     //calculate the mean free path
-     // by using normalized kinetic energy (= Ekin/mass)
-     G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
-     if ( rKineticEnergy > HighestBinValue) {
-       // beta >> 1
-       pathlength = ( rKineticEnergy + 1.0)* aCtau;
-     } else if ( rKineticEnergy > LowestBinValue) {
-       // check if aPhysicsTable exists
-       if (aPhysicsTable == 0) BuildPhysicsTable(*aParticleDef);
-       // beta is in the range valid for PhysicsTable
-       pathlength = aCtau *
-	 ((*aPhysicsTable)(0))-> GetValue(rKineticEnergy,isOutRange);
-     } else if ( rKineticEnergy < DBL_MIN ) {
-       // too slow particle
-#ifdef G4VERBOSE
-       if (GetVerboseLevel()>2) {
-	 G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
-         G4cout << aParticleDef->GetParticleName() << G4endl;
-	 G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
-       }
-#endif
-       pathlength = DBL_MIN;
-     } else {
-       // beta << 1
-       pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
-     }
-   }
-#ifdef G4VERBOSE
-   if (GetVerboseLevel()>1) {
-     G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
-   }
-#endif
-   return  pathlength;
+  } else {
+    //calculate the mean free path
+    // by using normalized kinetic energy (= Ekin/mass)
+    G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
+    if ( rKineticEnergy > HighestBinValue) {
+      // beta >> 1
+      pathlength = ( rKineticEnergy + 1.0)* aCtau;
+    } else if ( rKineticEnergy > LowestBinValue) {
+      // check if aPhysicsTable exists
+      if (aPhysicsTable == 0) BuildPhysicsTable(*aParticleDef);
+      // beta is in the range valid for PhysicsTable
+      pathlength = aCtau *
+	((*aPhysicsTable)(0))-> GetValue(rKineticEnergy,isOutRange);
+    } else if ( rKineticEnergy < DBL_MIN ) {
+      // too slow particle
+#ifdef G4VERBOSE
+      if (GetVerboseLevel()>2) {
+	G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
+	G4cout << aParticleDef->GetParticleName() << G4endl;
+	G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
+      }
+#endif
+      pathlength = DBL_MIN;
+    } else {
+      // beta << 1
+      pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
+    }
+  }
+#ifdef G4VERBOSE
+  if (GetVerboseLevel()>1) {
+    G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl;
+  }
+#endif
+  return  pathlength;
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -604,7 +609,7 @@
   G4int i;
   for ( i = 0 ; i < TotBin ; i++ ) {
-      gammainv = 1.0/(aVector->GetLowEdgeEnergy(i) + 1.0);
-      beta  = std::sqrt((1.0 - gammainv)*(1.0 +gammainv));
-      aVector->PutValue(i, beta/gammainv);
+    gammainv = 1.0/(aVector->GetLowEdgeEnergy(i) + 1.0);
+    beta  = std::sqrt((1.0 - gammainv)*(1.0 +gammainv));
+    aVector->PutValue(i, beta/gammainv);
   }
   aPhysicsTable->insert(aVector);
@@ -636,5 +641,5 @@
     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.");
+			      "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");
   }
   G4String dirName = getenv("G4RADIOACTIVEDATA");
@@ -647,5 +652,5 @@
   G4String file = os.str();
 
-  
+
   std::ifstream DecaySchemeFile(file);
   G4bool found(false);
@@ -664,5 +669,5 @@
       modeSumBR[i]       = 0.0;
     }
-    
+
     G4bool complete(false);
     char inputChars[80]={' '};
@@ -718,7 +723,7 @@
 	    a/=1000.;
 	    c/=1000.;
-		  
+
 	    //	cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl;
-		  
+
 	    switch (theDecayMode) {
 	    case IT:
@@ -750,5 +755,5 @@
 		    e0 = c*MeV/0.511;
 		    n = aBetaFermiFunction->GetFFN(e0);
-				
+
 		    // now to work out the histogram and initialise the random generator
 		    G4int npti = 100;				
@@ -765,5 +770,5 @@
 		      // Special treatment for K-40 (problem #1068) (flei,06/05/2010) 
 		      if (theParentNucleus.GetParticleName() == "K40[0.0]") f *=
-			(std::pow((g*g-1),3)+std::pow((ee-g),6)+7*(g*g-1)*std::pow((ee-g),2)*(g*g-1+std::pow((ee-g),2)));
+									      (std::pow((g*g-1),3)+std::pow((ee-g),6)+7*(g*g-1)*std::pow((ee-g),2)*(g*g-1+std::pow((ee-g),2)));
 		      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
 		    }		  
@@ -798,7 +803,7 @@
 		  if (e0 > 0.) {
 		    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, -(Z-1));
-				
+
 		    n = aBetaFermiFunction->GetFFN(e0);
-				
+
 		    // now to work out the histogram and initialise the random generator
 		    G4int npti = 100;				
@@ -824,5 +829,5 @@
 		    theDecayTable->Insert(aBetaPlusChannel);
 		    modeSumBR[2] += b;
-				
+
 		    delete[] pdf;
 		    delete aBetaFermiFunction;	      
@@ -909,5 +914,5 @@
 		G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
 			    FatalException, "Error in decay mode selection");
-		
+
 	      }
 	    }
@@ -972,6 +977,6 @@
 //
 void G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE, 
-				       G4int theG, std::vector<G4double> theRates, 
-				       std::vector<G4double> theTaos)
+				      G4int theG, std::vector<G4double> theRates, 
+				      std::vector<G4double> theTaos)
 { 
   //fill the decay rate vector 
@@ -993,5 +998,5 @@
   // 2) Add the coefficiencies to the decay rate table vector 
   //
-  
+
   //
   // Create and initialise variables used in the method.
@@ -999,10 +1004,11 @@
 
   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.);
   //
@@ -1012,7 +1018,8 @@
   G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy();
   G4double tao = theParentNucleus.GetPDGLifeTime();
+  if (tao < 0.) tao = 1e-30;
   taos.push_back(tao);
   G4int nEntry = 0;
-  
+
   //fill the decay rate with the intial isotope data
   SetDecayRate(Z,A,E,nGeneration,rates,taos);
@@ -1034,6 +1041,4 @@
   G4AlphaDecayChannel *theAlphaChannel = 0;
   G4RadioactiveDecayMode theDecayMode;
-  //  G4NuclearLevelManager levelManager;
-  //const G4NuclearLevel* level;
   G4double theBR = 0.0;
   G4int AP = 0;
@@ -1056,5 +1061,5 @@
   //
   theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
- 
+
   while (!stable) {
     nGeneration++;
@@ -1077,13 +1082,13 @@
 	aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus));
       }
+  	
+      G4DecayTable *theDecayTable = new G4DecayTable();
       aTempDecayTable = aParentNucleus->GetDecayTable();
+      for (i=0; i< 7; i++) brs[i] = 0.0;
+
       //
       //
       // Go through the decay table and to combine the same decay channels
       //
-      for (i=0; i< 7; i++) brs[i] = 0.0;
-      
-      G4DecayTable *theDecayTable = new G4DecayTable();
-      
       for (i=0; i<aTempDecayTable->entries(); i++) {
 	theChannel             = aTempDecayTable->GetDecayChannel(i);
@@ -1099,8 +1104,7 @@
 
 	  if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
-	    
-	    // Level hafe life is in ns and the gate as 1 micros by default
-	    if ( theDecayMode == 0 && level->HalfLife()*ns >= halflifethreshold ){    
-	      // some further though may needed here
+	    // 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);
 	    } 
@@ -1120,5 +1124,5 @@
       brs[3] = brs[4] =brs[5] =  0.0;
       for (i= 0; i<7; i++){
-	if (brs[i] > 0) {
+	if (brs[i] > 0.) {
 	  switch ( i ) {
 	  case 0:
@@ -1127,11 +1131,11 @@
 	    // Decay mode is isomeric transition.
 	    //
-	    
+
 	    theITChannel =  new G4ITDecayChannel
 	      (0, (const G4Ions*) aParentNucleus, brs[0]);
-	    
+
 	    theDecayTable->Insert(theITChannel);
 	    break;
-	    
+
 	  case 1:
 	    //
@@ -1142,7 +1146,7 @@
 							       brs[1], 0.*MeV, 0.*MeV, 1, false, 0);
 	    theDecayTable->Insert(theBetaMinusChannel);
-	    
+
 	    break;
-	    
+
 	  case 2:
 	    //
@@ -1154,5 +1158,5 @@
 	    theDecayTable->Insert(theBetaPlusChannel);
 	    break;		      
-	    
+
 	  case 6:
 	    //
@@ -1164,5 +1168,5 @@
 	    theDecayTable->Insert(theAlphaChannel);
 	    break;
-	    
+
 	  default:
 	    break;
@@ -1170,5 +1174,4 @@
 	}
       }
-	
       // 
       // loop over all braches in theDecayTable
@@ -1179,10 +1182,11 @@
 	theBR = theChannel->GetBR();
 	theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus();
-	
-	//
-	// now test if the daughterNucleus is a valid one
-	//
-	if (IsApplicable(*theDaughterNucleus) && theBR 
-	    && aParentNucleus != theDaughterNucleus ) {
+	//  first check if the decay of the original nucleus is an IT channel, if true create a new groud-level nucleus
+	if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1 ) {
+	    A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass();
+	    Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
+	    theDaughterNucleus=theIonTable->GetIon(Z,A,0.);
+	}
+	if (IsApplicable(*theDaughterNucleus) && theBR && aParentNucleus != theDaughterNucleus) { 
 	  // need to make sure daugher has decaytable
 	  if (!IsLoaded(*theDaughterNucleus)){
@@ -1194,32 +1198,45 @@
 	    Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber();
 	    E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy();
-	  
+
 	    TaoPlus = theDaughterNucleus->GetPDGLifeTime();
 	    //		cout << TaoPlus <<G4endl;
-	    if (TaoPlus > 0.) {
-	      // first set the taos, one simply need to add to the parent ones
-	      taos.clear();
-	      taos = TP;
-	      taos.push_back(TaoPlus);
-	      // now calculate the coefficiencies
-	      //
-	      // they are in two parts, first the les than n ones
-	      rates.clear();
-	      size_t k;
-	      for (k = 0; k < RP.size(); k++){
+	    if (TaoPlus <= 0.)  TaoPlus = 1e-30; 
+
+	    // first set the taos, one simply need to add to the parent ones
+	    taos.clear();
+	    taos = TP;
+	    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();
+	    size_t k;
+	    for (k = 0; k < RP.size(); k++){
+	      if ((TP[k]-TaoPlus) == 0.) {
+		theRate =  1e30;
+	      } else {
 		theRate = TP[k]/(TP[k]-TaoPlus) * theBR * RP[k];
-		rates.push_back(theRate);
 	      }
-	      //
-	      // the sencond part: the n:n coefficiency
-	      theRate = 0.;
-	      for (k = 0; k < RP.size(); k++){
-		theRate -=TaoPlus/(TP[k]-TaoPlus) * 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.;
+	    G4double aRate;
+	    for (k = 0; k < RP.size(); k++){
+	      if ((TP[k]-TaoPlus) == 0.) {
+		aRate =  1e30;
+	      } else {
+		aRate = TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k];
 	      }
-	      rates.push_back(theRate); 	      
-	      SetDecayRate (Z,A,E,nGeneration,rates,taos);
-	      theDecayRateVector.push_back(theDecayRate);
-	      nEntry++;
-	    } 
+	      theRate -= aRate;
+	    }
+	    rates.push_back(theRate); 	      
+	    SetDecayRate (Z,A,E,nGeneration,rates,taos);
+	    theDecayRateVector.push_back(theDecayRate);
+	    nEntry++;
 	  }
 	}  
@@ -1233,9 +1250,9 @@
     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) 
@@ -1246,7 +1263,7 @@
   // now fill the decay table with the newly completed decay rate vector
   //
-  
+
   theDecayRateTable.SetItsRates(theDecayRateVector);
-  
+
   //
   // finally add the decayratetable to the tablevector
@@ -1254,5 +1271,5 @@
   theDecayRateTableVector.push_back(theDecayRateTable);
 }
-  
+
 ////////////////////////////////////////////////////////////////////////////////
 //
@@ -1267,6 +1284,6 @@
   std::ifstream infile ( filename, std::ios::in );
   if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open source data file" );
-  
-  float bin, flux;
+
+  G4double bin, flux;
   NSourceBin = -1;
   while (infile >> bin >> flux ) {
@@ -1277,4 +1294,6 @@
   }
   SetAnalogueMonteCarlo(0);
+  infile.close();
+
 #ifdef G4VERBOSE
   if (GetVerboseLevel()>1)
@@ -1294,6 +1313,12 @@
   std::ifstream infile ( filename, std::ios::in);
   if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException,  "Unable to open bias data file" );
-  
-  float bin, flux;
+
+  G4double bin, flux;
+  G4int dWindows = 0;
+  G4int i ;
+
+  theRadioactivityTables.clear();
+  //  for (i = 0; i<100; i++) decayWindows[i] = -1;
+
   NDecayBin = -1;
   while (infile >> bin >> flux ) {
@@ -1302,9 +1327,18 @@
     DBin[NDecayBin] = bin * s;
     DProfile[NDecayBin] = flux;
+    if (flux > 0.) {
+      decayWindows[NDecayBin] = dWindows;
+      dWindows++;
+      G4RadioactivityTable *rTable = new G4RadioactivityTable() ;
+      theRadioactivityTables.push_back(rTable);
+    }
   }
-  G4int i ;
   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)
@@ -1315,13 +1349,12 @@
 ////////////////////////////////////////////////////////////////////////////////
 //
-//
 // DecayIt
 //
-G4VParticleChange* G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step& )
-{
-  //
+G4VParticleChange*
+G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&)
+{
   // Initialize the G4ParticleChange object. Get the particle details and the
   // decay table.
-  //
+
   fParticleChangeForRadDecay.Initialize(theTrack);
   const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle();
@@ -1329,58 +1362,51 @@
 
   // First check whether RDM applies to the current logical volume
-  //
-  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;
+
+  if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(), 
+			  theTrack.GetVolume()->GetLogicalVolume()->GetName())) {
+#ifdef G4VERBOSE
+    if (GetVerboseLevel()>0) {
+      G4cout <<"G4RadioactiveDecay::DecayIt : "
+             << theTrack.GetVolume()->GetLogicalVolume()->GetName()
+             << " is not selected for the RDM"<< G4endl;
+      G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl;
+      G4cout << " The Valid volumes are " << G4endl;
+      for (size_t i = 0; i< ValidVolumes.size(); i++) G4cout << ValidVolumes[i] << G4endl;
     }
-   
+#endif
+    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
+
+    // Kill the parent particle.
+
+    fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
+    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
+    ClearNumberOfInteractionLengthLeft();
+    return &fParticleChangeForRadDecay;
+  }
+
   // now check is the particle is valid for RDM
-  //
-  if (!(IsApplicable(*theParticleDef)))
-    { 
-      //
-      // The particle is not a Ion or outside the nucleuslimits for decay
-      //
-#ifdef G4VERBOSE
-      if (GetVerboseLevel()>0)
-	{
-	  G4cerr <<"G4RadioactiveDecay::DecayIt : "
-		 <<theParticleDef->GetParticleName() 
-		 << " is not a valid nucleus for the RDM"<< G4endl;
-	}
-#endif
-      fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
-      //
-      //
-      // Kill the parent particle.
-      //
-      fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
-      fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
-      ClearNumberOfInteractionLengthLeft();
-      return &fParticleChangeForRadDecay;
+
+  if (!(IsApplicable(*theParticleDef))) { 
+    //
+    // The particle is not a Ion or outside the nucleuslimits for decay
+    //
+#ifdef G4VERBOSE
+    if (GetVerboseLevel()>0) {
+      G4cerr <<"G4RadioactiveDecay::DecayIt : "
+	     <<theParticleDef->GetParticleName() 
+	     << " is not a valid nucleus for the RDM"<< G4endl;
     }
-  
+#endif
+    fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
+
+    //
+    // Kill the parent particle.
+    //
+    fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
+    fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
+    ClearNumberOfInteractionLengthLeft();
+    return &fParticleChangeForRadDecay;
+  }
+
   if (!IsLoaded(*theParticleDef))
     {
@@ -1388,5 +1414,5 @@
     }
   G4DecayTable *theDecayTable = theParticleDef->GetDecayTable();
-  
+
   if  (theDecayTable == 0 || theDecayTable->entries() == 0 )
     {
@@ -1423,5 +1449,5 @@
       G4ThreeVector currentPosition;
       currentPosition = theTrack.GetPosition();
-      
+
       // check whether use Analogue or VR implementation
       //
@@ -1454,5 +1480,5 @@
 	G4double ParentEnergy = theParticle->GetTotalEnergy();
 	G4ThreeVector ParentDirection(theParticle->GetMomentumDirection());
-	
+
 	if (theTrack.GetTrackStatus() == fStopButAlive )
 	  {
@@ -1501,5 +1527,5 @@
 	      (products->PopProducts(), finalGlobalTime, currentPosition);
 	    secondary->SetGoodForTrackingFlag();
-			secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
+	    secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
 	    fParticleChangeForRadDecay.AddSecondary(secondary);
 	  }
@@ -1521,5 +1547,5 @@
 	if (!IsRateTableReady(*theParticleDef)) {
 	  // if the decayrates are not ready, calculate them and 
-  	  // add to the rate table vector 
+	  // add to the rate table vector 
 	  AddDecayRateTable(*theParticleDef);
 	}
@@ -1538,5 +1564,5 @@
 	G4double taotime;
 	G4double decayRate;
-	
+
 	size_t i;
 	size_t j;
@@ -1557,30 +1583,28 @@
 	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);
 	    
-	    // claculate the first part of the weight function
+	    // calculate the first part of the weight function
 	    
-	    G4double weight1 =1./DProfile[nbin-1] 
-	      *(DBin[nbin]-DBin[nbin-1])
-	      /NSplit;
-	    if (nbin > 1) {
-	       weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
-		 *(DBin[nbin]-DBin[nbin-1])
-		 /NSplit;}
+	    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++){
+
+	    for (i = 0; i<theDecayRateVector.size(); i++){
 	      PZ = theDecayRateVector[i].GetZ();
 	      PA = theDecayRateVector[i].GetA();
@@ -1589,53 +1613,42 @@
 	      PR = theDecayRateVector[i].GetDecayRateC();
 
-	      //
-	      // Get the decay time following the decay probability function 
-	      // suppllied by user
-	      //
-	      G4double theDecayTime = GetDecayTime();
-	      
-	      G4int nbin = GetDecayTimeBin(theDecayTime);
-	      
-	      // claculate the first part of the weight function
-	      
-	      G4double weight1 =1./DProfile[nbin-1] 
-		*(DBin[nbin]-DBin[nbin-1])
-		/NSplit;
-	      if (nbin > 1) {
-		weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])
-		  *(DBin[nbin]-DBin[nbin-1])
-		  /NSplit;}
-	      // it should be calculated in seconds
-	      weight1 /= s ;
-	      
 	      // a temprary products buffer and its contents is transfered to 
 	      // the products at the end of the loop
-	      //
+
 	      G4DecayProducts *tempprods = 0;
-	      
+
 	      // calculate the decay rate of the isotope
-	      // one need to fold the the source bias function with the decaytime
-	      //
+	      // 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.;
 	      for ( j = 0; j < PT.size(); j++){
 		taotime = GetTaoTime(theDecayTime,PT[j]);
 		decayRate -= PR[j] * taotime;
+		// Eq.4.23 of of the TN
+		// note the negative here is required as the rate in the eqation is defined to be negative, 
+		// i.e. decay away, but we need pasitive value here.
+
+		//		G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t" << decayRate << G4endl ;		
 	      }
-	      
-	      // decayRatehe 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
-	      
-	      
-	      //
+
 	      // now calculate the statistical weight
-	      //
-	      
-	      G4double weight = weight1*decayRate; 
+
+	      // one need to fold the 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(); 
+		
+	      // 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,weight);
+				
 	      // decay the isotope 
 	      theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable());
 	      parentNucleus = theIonTable->GetIon(PZ,PA,PE);
-	      
+
 	      // decide whther to apply branching ratio bias or not
 	      //
@@ -1644,49 +1657,44 @@
 		ndecaych = G4int(theDecayTable->entries()*G4UniformRand());
 		G4VDecayChannel *theDecayChannel = theDecayTable->GetDecayChannel(ndecaych);
-		if (theDecayChannel == 0)
-		  {
-		    // Decay channel not found.
-#ifdef G4VERBOSE
-		    if (GetVerboseLevel()>0)
-		      {
-			G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
-			G4cerr <<G4endl;
-			theDecayTable ->DumpInfo();
-		      }
-#endif
+		if (theDecayChannel == 0) {
+		  // Decay channel not found.
+#ifdef G4VERBOSE
+		  if (GetVerboseLevel()>0) {
+		    G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
+		    G4cerr <<G4endl;
+		    theDecayTable ->DumpInfo();
 		  }
-		else
-		  {
-		    // A decay channel has been identified, so execute the DecayIt.
-		    G4double tempmass = parentNucleus->GetPDGMass();      
-		    tempprods = theDecayChannel->DecayIt(tempmass);
-		    weight *= (theDecayChannel->GetBR())*(theDecayTable->entries());
-		  }
-	      }
-	      else {
+#endif
+		} else {
+		  // A decay channel has been identified, so execute the DecayIt.
+		  G4double tempmass = parentNucleus->GetPDGMass();      
+		  tempprods = theDecayChannel->DecayIt(tempmass);
+		  weight *= (theDecayChannel->GetBR())*(theDecayTable->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);
-		  }
+	      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
@@ -1699,7 +1707,7 @@
 					     secondaryparticles[index], ptime[index], currentPosition);
 	    secondary->SetGoodForTrackingFlag(); 	   
-			secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
+	    secondary->SetTouchableHandle(theTrack.GetTouchableHandle());
 	    secondary->SetWeight(pw[index]); 	   
-            fParticleChangeForRadDecay.AddSecondary(secondary); 
+	    fParticleChangeForRadDecay.AddSecondary(secondary); 
 	  }
 	//
@@ -1708,7 +1716,7 @@
 	//theTrack.SetTrackStatus(fStopButAlive);
 	//energyDeposit += theParticle->GetKineticEnergy();
-	
+
       }
-    
+
       //
       // Kill the parent particle.
@@ -1722,76 +1730,54 @@
       //
       ClearNumberOfInteractionLengthLeft();
-      
+
       return &fParticleChangeForRadDecay ;
     }
 } 
 
-////////////////////////////////////////////////////////////////////////////////
-//
+///////////////////////////////////////////////////////////////////
 //
 // DoDecay
 //
-G4DecayProducts* G4RadioactiveDecay::DoDecay(  G4ParticleDefinition& theParticleDef )
-{
-  G4DecayProducts *products = 0;
-  //
-  //
+G4DecayProducts*
+G4RadioactiveDecay::DoDecay(  G4ParticleDefinition& theParticleDef )
+{
+  G4DecayProducts* products = 0;
+
   // follow the decaytable and generate the secondaries...
-  // 
-  //
-#ifdef G4VERBOSE
-  if (GetVerboseLevel()>0)
-    {
-      G4cout<<"Begin of DoDecay..."<<G4endl;
+
+#ifdef G4VERBOSE
+  if (GetVerboseLevel()>0) G4cout<<"Begin of DoDecay..."<<G4endl;
+#endif
+
+  G4DecayTable* theDecayTable = theParticleDef.GetDecayTable();
+
+  // Choose a decay channel.
+
+#ifdef G4VERBOSE
+  if (GetVerboseLevel()>0) G4cout <<"Selecte a channel..."<<G4endl;
+#endif
+
+  G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel();
+  if (theDecayChannel == 0) {
+    // Decay channel not found.
+
+    G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
+    G4cerr <<G4endl;
+    theDecayTable ->DumpInfo();
+  } else {
+
+    // A decay channel has been identified, so execute the DecayIt.
+
+#ifdef G4VERBOSE
+    if (GetVerboseLevel()>1) {
+      G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
+      G4cerr <<theDecayChannel <<G4endl;
     }
 #endif
-  G4DecayTable *theDecayTable = theParticleDef.GetDecayTable();
-  //
-  // Choose a decay channel.
-  //
-#ifdef G4VERBOSE
-  if (GetVerboseLevel()>0)
-    {
-      G4cout <<"Selecte a channel..."<<G4endl;
-    }
-#endif
-  G4VDecayChannel *theDecayChannel = theDecayTable->SelectADecayChannel();
-  if (theDecayChannel == 0)
-    {
-      // Decay channel not found.
-      //
-      G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel";
-      G4cerr <<G4endl;
-      theDecayTable ->DumpInfo();
-    }
-      else
-    {
-      //
-      // A decay channel has been identified, so execute the DecayIt.
-      //
-#ifdef G4VERBOSE
-      if (GetVerboseLevel()>1)
-	{
-	  G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel  addr:";
-	  G4cerr <<theDecayChannel <<G4endl;
-	}
-#endif
-      
-      G4double tempmass = theParticleDef.GetPDGMass();
-      //
-      
-      products = theDecayChannel->DecayIt(tempmass);
-      
-    }
+
+    G4double tempmass = theParticleDef.GetPDGMass();
+    products = theDecayChannel->DecayIt(tempmass);
+  }
+
   return products;
-
-}
-
-
-
-
-
-
-
-
-
+}
Index: trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecaymessenger.cc
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecaymessenger.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecaymessenger.cc	(revision 1347)
@@ -123,5 +123,5 @@
   icmCmd->SetParameterName("applyICM",true);
   icmCmd->SetDefaultValue(true);
-  icmCmd->AvailableForStates(G4State_PreInit);
+  //icmCmd->AvailableForStates(G4State_PreInit);
   //
   // Command contols whether ARM will be applied or not
@@ -131,5 +131,5 @@
   armCmd->SetParameterName("applyARM",true);
   armCmd->SetDefaultValue(true);
-  armCmd->AvailableForStates(G4State_PreInit);
+  //armCmd->AvailableForStates(G4State_PreInit);
   //
   // Command to set the h-l thresold for isomer production
@@ -140,5 +140,5 @@
   // hlthCmd->SetRange("hlThreshold>0.");
   hlthCmd->SetUnitCategory("Time");
-  hlthCmd->AvailableForStates(G4State_PreInit);
+  //  hlthCmd->AvailableForStates(G4State_PreInit);
   //
   // Command to define the incident particle source time profile.
@@ -209,15 +209,7 @@
 				    SetNucleusLimits(nucleuslimitsCmd->GetNewNucleusLimitsValue(newValues));}
   else if  (command==analoguemcCmd) {
-    G4int vl;
-    const char* t = newValues;
-    std::istringstream is(t);
-    is >> vl;
-    theRadioactiveDecayContainer->SetAnalogueMonteCarlo(vl!=0);}
+    theRadioactiveDecayContainer->SetAnalogueMonteCarlo(analoguemcCmd->GetNewBoolValue(newValues));}
   else if  (command==fbetaCmd) {
-    G4int vl;
-    const char* t = newValues;
-    std::istringstream is(t);
-    is >> vl;
-    theRadioactiveDecayContainer->SetFBeta(vl!=0);}
+    theRadioactiveDecayContainer->SetFBeta(fbetaCmd->GetNewBoolValue(newValues));}
   else if  (command==avolumeCmd) {theRadioactiveDecayContainer->
 				   SelectAVolume(newValues);}
@@ -229,9 +221,5 @@
 				   DeselectAllVolumes();}
   else if  (command==brbiasCmd) {
-    G4int vl;
-    const char* t = newValues;
-    std::istringstream is(t);
-    is >> vl;
-    theRadioactiveDecayContainer->SetBRBias(vl!=0);}
+    theRadioactiveDecayContainer->SetBRBias(brbiasCmd->GetNewBoolValue(newValues));}
   else if (command==sourcetimeprofileCmd) {theRadioactiveDecayContainer->
 					     SetSourceTimeProfile(newValues);}
