Index: trunk/source/processes/hadronic/models/radioactive_decay/History
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/History	(revision 1201)
+++ trunk/source/processes/hadronic/models/radioactive_decay/History	(revision 1315)
@@ -15,8 +15,38 @@
      ----------------------------------------------------------
 
+26 May 2010 Dennis Wright (radioactive_decay-V09-03-00)
+-------------------------------------------------------
+tag HEAD to submit for testing
+
+11 may 2010 F.Lei
+- G4RadioactiveDecay.cc::LoadDecayTable:line 931 create a decaytable for isomers not included in RDM database and assume 
+  they all under go IT decay.
+- G4RadioactiveDecay.cc: line 1434 after DoDecay() check if there is any secondary produced. Kill the track if not to 
+  prevent it entering a infinite loop
+- G4RadioactiveDecay.cc: line 1459. make sure the propertime is positive. negative case occures when the isomer is not in
+  RDM database and its f-l is set to -1 by default.  
+- G4NuclearDecayChannel.cc: Added initialisation to all three constructors: 
+  halflifethreshold = 1e-6*second;
+  applyICM = true;
+  applyARM = true;
+
+10 May 2010 F.Lei
+- G4RadioactiveDecay.hh:line 264 corrected typo as pointed out by Luciano Pandola
+- G4RadioactiveDecay.cc: insert a special treatment for K-40 beta decay at line 774 as suggested by Mauro Taiuti.
+  This fixes bug #1068
+  applied SetICM() SetARM() and SetHLThreshold() to a decay channel is created.
+- G4NuclearDecayChannel.cc: line 396. Limit the shell index to < 7, as required by the ARM.
+  line 329: change to use BreakUp rather than BreakItUp so to limit to one transition per step when ICM is applied (bug #1001)
+  (Note: changes have also been made in 30/08/2008 to use the correct atomic shell index and to conserve energy atfer ARM, 
+   these problems were pointed out by Andreas Zoglauer in the bug report).
+- G4RadioactiveDecay.hh & G4NuclearDecayChannel.hh: added private member data "applyICM" "applyARM" "halflifethresold" 
+  and their public Set methods SetICM(), SetARM() and SetHLThreshold().
+- G4RadioactiveDecaymessenger: added the UICommands for SetICM, SETARM and SetHLThreshold.
+
 24 July 2009 V.Ivanchenko (radioactive_decay-V09-02-00)
 - GNUmakefile - added dependence on electromagnetic/utils
 
 09 July 2008 Dennis Wright (radioactive_decay-V09-01-02)
+--------------------------------------------------------
 - replace exit() with G4Exception in G4RadioactiveDecay and G4NuclearDecayChannel
 
Index: trunk/source/processes/hadronic/models/radioactive_decay/include/G4NuclearDecayChannel.hh
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/include/G4NuclearDecayChannel.hh	(revision 1201)
+++ trunk/source/processes/hadronic/models/radioactive_decay/include/G4NuclearDecayChannel.hh	(revision 1315)
@@ -88,4 +88,13 @@
   // Returns the decay products
   //
+  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
+  // 
   inline G4RadioactiveDecayMode GetDecayMode () {return decayMode;}
   // Returns the decay mode
@@ -136,4 +145,7 @@
   G4double               FermiFN;
   G4bool                 BetaSimple;
+  G4double               halflifethreshold;
+  G4bool                 applyICM;
+  G4bool                 applyARM;
   CLHEP::RandGeneral*    RandomEnergy;    
 };
Index: trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecay.hh
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecay.hh	(revision 1201)
+++ trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecay.hh	(revision 1315)
@@ -114,4 +114,13 @@
   //   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"
@@ -235,4 +244,8 @@
   G4int                         NSplit;
 
+  G4double                      halflifethreshold;
+  G4bool                        applyICM;
+  G4bool                        applyARM;
+
   G4int                         NSourceBin;
   G4double                      SBin[99];
@@ -243,6 +256,6 @@
   G4double                      DProfile[99];
 
-  std::vector<G4String>              LoadedNuclei;
-  std::vector<G4String>              ValidVolumes;
+  std::vector<G4String>         LoadedNuclei;
+  std::vector<G4String>         ValidVolumes;
 
   G4RadioactiveDecayRate        theDecayRate;
@@ -262,5 +275,5 @@
   G4ParticleChangeForRadDecay   fParticleChangeForRadDecay;
 
-  inline G4double AtRestGetPhysicalInteractionLengt
+  inline G4double AtRestGetPhysicalInteractionLength
     (const G4Track& track, G4ForceCondition* condition)
   {fRemainderLifeTime = G4VRestDiscreteProcess::
Index: trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecaymessenger.hh
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecaymessenger.hh	(revision 1201)
+++ trunk/source/processes/hadronic/models/radioactive_decay/include/G4RadioactiveDecaymessenger.hh	(revision 1315)
@@ -107,4 +107,7 @@
   G4UIcmdWithoutParameter        *allvolumesCmd;
   G4UIcmdWithoutParameter        *deallvolumesCmd;
+  G4UIcmdWithABool               *icmCmd;
+  G4UIcmdWithABool               *armCmd;
+  G4UIcmdWithADoubleAndUnit      *hlthCmd;
 };
 
Index: trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc	(revision 1201)
+++ trunk/source/processes/hadronic/models/radioactive_decay/src/G4NuclearDecayChannel.cc	(revision 1315)
@@ -101,4 +101,7 @@
   FillDaughterNucleus (0, A, Z, theDaughterExcitation);
   Qtransition = theQtransition;
+  halflifethreshold = 1e-6*second;
+  applyICM = true;
+  applyARM = true;
 }
 ///////////////////////////////////////////////////////////////////////////////
@@ -131,4 +134,7 @@
   FillDaughterNucleus (1, A, Z, theDaughterExcitation);
   Qtransition = theQtransition;
+  halflifethreshold = 1e-6*second;
+  applyICM = true;
+  applyARM = true;
 }
 ///////////////////////////////////////////////////////////////////////////////
@@ -171,4 +177,7 @@
   Qtransition = theQtransition;
   FermiFN = theFFN;
+  halflifethreshold = 1e-6*second;
+  applyICM = true;
+  applyARM = true;
 }
 
@@ -288,10 +297,10 @@
   //
   // If the decay is to an excited state of the daughter nuclide, we need
-  // to apply the photo-evaporation process.
+  // to apply the photo-evaporation process. This includes the IT decay mode itself.
   //
   // needed to hold the shell idex after ICM
   G4int shellIndex = -1;
   //
-  if (daughterExcitation > 0.0)
+  if (daughterExcitation > 0.0) 
     {
       //
@@ -311,9 +320,9 @@
       deexcitation->SetVerboseLevel(GetVerboseLevel());
       // switch on/off internal electron conversion
-      deexcitation->SetICM(true);
+      deexcitation->SetICM(applyICM);
       // set the maximum life-time for a level that will be treated. Level with life-time longer than this
       // will be outputed as meta-stable isotope
       //
-      deexcitation->SetMaxHalfLife(1e-6*second);
+      deexcitation->SetMaxHalfLife(halflifethreshold);
       // but in IT mode, we need to force the transition 
       if (decayMode == 0) {
@@ -323,8 +332,14 @@
       }
       //
-      // Get the gammas by deexciting the nucleus.
-      //
-      G4FragmentVector* gammas = deexcitation->BreakItUp(nucleus);
-      // in the case of BreakItUp(nucleus), the returned G4FragmentVector contains the residual nuclide
+      // Now apply the photo-evaporation
+      // Use BreakUp() so limit to one transition at a time, if ICM is requested
+      // this change is realted to bug#1001  (F.Lei 07/05/2010)
+      G4FragmentVector* gammas = 0;	
+      if (applyICM) {
+	gammas = deexcitation->BreakUp(nucleus);	
+      } else {
+	gammas = deexcitation->BreakItUp(nucleus);
+      }
+      // the returned G4FragmentVector contains the residual nuclide
       // as its last entry.
       G4int nGammas=gammas->size()-1;
@@ -343,5 +358,5 @@
 	}
       //
-      //      now the nucleus
+      // now the nucleus
       G4double finalDaughterExcitation = gammas->operator[](nGammas)->GetExcitationEnergy();
       // f.lei (03/01/03) this is needed to fix the crach in test18 
@@ -356,4 +371,5 @@
 	 daughterMomentum1);
       products->PushProducts (dynamicDaughter); 
+      
       // retrive the ICM shell index
       shellIndex = deexcitation->GetVacantShellNumber();
@@ -392,5 +408,7 @@
 	//
 	{
-	  eShell = G4int(G4UniformRand()*5)+4;
+	// limit the shell index to 6 as specified by the ARM (F.Lei 06/05/2010)
+	// eShell = G4int(G4UniformRand()*5)+4;
+	  eShell = G4int(G4UniformRand()*3)+4;
 	}
 	break;
@@ -401,12 +419,13 @@
       }
   }
-  // now deal with the IT case where ICM may have been applied
+  // Also have to deal with the IT case where ICM may have been applied
   //
   if (decayMode == 0) {
     eShell = shellIndex;
   }
-  // now apply ARM if there is a vaccancy
-  //
-  if (eShell != -1) {
+  //
+  // now apply ARM if it is requested and there is a vaccancy
+  //
+  if (applyARM && eShell != -1) {
     G4int aZ = daughterZ;
     if (aZ > 5 && aZ < 100) {  // only applies to 5< Z <100 
Index: trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc	(revision 1201)
+++ trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc	(revision 1315)
@@ -168,4 +168,7 @@
   FBeta       = false ;
   BRBias      = true ;
+  applyICM    = true ;
+  applyARM    = true ;
+  halflifethreshold = 1e-6*second;
   //
   // RDM applies to xall logical volumes as default
@@ -200,7 +203,7 @@
   // All particles, other than G4Ions, are rejected by default.
   //
+  if (((const G4Ions*)(&aParticle))->GetExcitationEnergy() > 0.) {return true;}
   if (aParticle.GetParticleName() == "GenericIon")      {return true;}
   else if (!(aParticle.GetParticleType() == "nucleus") || aParticle.GetPDGLifeTime() < 0. ) {return false;}
-
   //
   //
@@ -482,4 +485,6 @@
     else if (theLife < 0.0) {meanlife = DBL_MAX;}
     else {meanlife = theLife;}
+  // set the meanlife to zero for excited istopes which is not in the RDM database
+    if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;}
   }
 #ifdef G4VERBOSE
@@ -612,6 +617,5 @@
 //     theParentNucleus.
 //
-G4DecayTable *G4RadioactiveDecay::LoadDecayTable (G4ParticleDefinition
-  &theParentNucleus)
+G4DecayTable *G4RadioactiveDecay::LoadDecayTable (G4ParticleDefinition &theParentNucleus)
 {
   //
@@ -645,303 +649,304 @@
   
   std::ifstream DecaySchemeFile(file);
-
-  if (!DecaySchemeFile)
+  G4bool found(false);
+  if (DecaySchemeFile) { 
     //
+    //
+    // Initialise variables used for reading in radioactive decay data.
+    //
+    G4int    nMode = 7;
+    G4bool   modeFirstRecord[7];
+    G4double modeTotalBR[7];
+    G4double modeSumBR[7];
+    G4int i;
+    for (i=0; i<nMode; i++) {
+      modeFirstRecord[i] = true;
+      modeSumBR[i]       = 0.0;
+    }
+    
+    G4bool complete(false);
+    char inputChars[80]={' '};
+    G4String inputLine;
+    G4String recordType("");
+    G4RadioactiveDecayMode theDecayMode;
+    G4double a(0.0);
+    G4double b(0.0);
+    G4double c(0.0);
+    G4double n(1.0);
+    G4double e0;
+    //
+    //
+    // Go through each record in the data file until you identify the decay
+    // data relating to the nuclide of concern.
+    //
+    //  while (!complete && -DecaySchemeFile.getline(inputChars, 80).eof() != EOF)
+    while (!complete && !DecaySchemeFile.getline(inputChars, 80).eof()) {
+      inputLine = inputChars;
+      //    G4String::stripType stripend(1);
+      //    inputLine = inputLine.strip(stripend);
+      inputLine = inputLine.strip(1);
+      if (inputChars[0] != '#' && inputLine.length() != 0) {
+	std::istringstream tmpStream(inputLine);
+	if (inputChars[0] == 'P') {
+	  //
+	  //
+	  // Nucleus is a parent type.  Check the excitation level to see if it matches
+	  // that of theParentNucleus
+	  //
+	  tmpStream >>recordType >>a >>b;
+	  if (found) {complete = true;}
+	  else {found = (std::abs(a*keV - E)<levelTolerance);}
+	} else if (found) {
+	  //
+	  //
+	  // The right part of the radioactive decay data file has been found.  Search
+	  // through it to determine the mode of decay of the subsequent records.
+	  //
+	  if (inputChars[0] == 'W') {
+#ifdef G4VERBOSE
+	    if (GetVerboseLevel()>0) {
+	      // a comment line identified and print out the message
+	      //
+	      G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
+	      G4cout << "   In data file " << file << G4endl;
+	      G4cout << "   " << inputLine << G4endl;
+	    }
+#endif
+	  }	
+	  else {
+	    tmpStream >>theDecayMode >>a >>b >>c;
+	    a/=1000.;
+	    c/=1000.;
+		  
+	    //	cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl;
+		  
+	    switch (theDecayMode) {
+	    case IT:
+	      {
+		//
+		//
+		// Decay mode is isomeric transition.
+		//
+		G4ITDecayChannel *anITChannel = new G4ITDecayChannel
+		  (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, b);
+		anITChannel->SetICM(applyICM);
+		anITChannel->SetARM(applyARM);
+		anITChannel->SetHLThreshold(halflifethreshold);
+		theDecayTable->Insert(anITChannel);
+		break;
+	      }
+	    case BetaMinus:
+	      {
+		//
+		//
+		// Decay mode is beta-.
+		//
+		if (modeFirstRecord[1])
+		  {modeFirstRecord[1] = false; modeTotalBR[1] = b;}
+		else {
+		  if (c > 0.) {
+		    // to work out the Fermi function normalization factor first
+		    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, (Z+1));
+		    e0 = c*MeV/0.511;
+		    n = aBetaFermiFunction->GetFFN(e0);
+				
+		    // now to work out the histogram and initialise the random generator
+		    G4int npti = 100;				
+		    G4double* pdf = new G4double[npti];
+		    G4int ptn;
+		    G4double g,e,ee,f;
+		    ee = e0+1.;
+		    for (ptn=0; ptn<npti; ptn++) {
+		      // e =e0*(ptn+1.)/102.;
+		      // bug fix (#662) (flei, 22/09/2004)
+		      e =e0*(ptn+0.5)/100.;
+		      g = e+1.;
+		      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
+		      // 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)));
+		      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
+		    }		  
+		    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
+		    G4BetaMinusDecayChannel *aBetaMinusChannel = new
+		      G4BetaMinusDecayChannel (GetVerboseLevel(), &theParentNucleus,
+					       b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
+		    theDecayTable->Insert(aBetaMinusChannel);
+		    modeSumBR[1] += b;
+		    delete[] pdf;
+		    delete aBetaFermiFunction;
+		  }
+		}
+		break;
+	      case BetaPlus:
+		//
+		//
+		// Decay mode is beta+.
+		//
+		if (modeFirstRecord[2])
+		  {modeFirstRecord[2] = false; modeTotalBR[2] = b;}
+		else {
+		  //  e0 = c*MeV/0.511;
+		  // bug fix (#662) (flei, 22/09/2004)
+		  // need to test e0 as there are some data files (e.g. z67.a162) which have entries for beta+
+		  // with Q < 2Me
+		  //
+		  e0 = c*MeV/0.511 -2.;
+		  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;				
+		    G4double* pdf = new G4double[npti];
+		    G4int ptn;
+		    G4double g,e,ee,f;
+		    ee = e0+1.;
+		    for (ptn=0; ptn<npti; ptn++) {
+		      // e =e0*(ptn+1.)/102.;
+		      // bug fix (#662) (flei, 22/09/2004)
+		      e =e0*(ptn+0.5)/100.;
+		      g = e+1.;
+		      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
+		      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
+		    }
+		    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
+		    G4BetaPlusDecayChannel *aBetaPlusChannel = new 
+		      G4BetaPlusDecayChannel (GetVerboseLevel(), &theParentNucleus,
+					      b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
+		    theDecayTable->Insert(aBetaPlusChannel);
+		    modeSumBR[2] += b;
+				
+		    delete[] pdf;
+		    delete aBetaFermiFunction;	      
+		  }
+		}
+		break;
+	      case KshellEC:
+		//
+		//
+		// Decay mode is K-electron capture.
+		//
+		if (modeFirstRecord[3])
+		  {modeFirstRecord[3] = false; modeTotalBR[3] = b;}
+		else {
+		  G4KshellECDecayChannel *aKECChannel = new
+		    G4KshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
+					    b, c*MeV, a*MeV);
+		  aKECChannel->SetICM(applyICM);
+		  aKECChannel->SetARM(applyARM);
+		  aKECChannel->SetHLThreshold(halflifethreshold);
+		  theDecayTable->Insert(aKECChannel);
+		  //delete aKECChannel;
+		  modeSumBR[3] += b;
+		}
+		break;
+	      case LshellEC:
+		//
+		//
+		// Decay mode is L-electron capture.
+		//
+		if (modeFirstRecord[4])
+		  {modeFirstRecord[4] = false; modeTotalBR[4] = b;}
+		else {
+		  G4LshellECDecayChannel *aLECChannel = new
+		    G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
+					    b, c*MeV, a*MeV);
+		  aLECChannel->SetICM(applyICM);
+		  aLECChannel->SetARM(applyARM);
+		  aLECChannel->SetHLThreshold(halflifethreshold);
+		  theDecayTable->Insert(aLECChannel);
+		  //delete aLECChannel;
+		  modeSumBR[4] += b;
+		}
+		break;
+	      case MshellEC:
+		//
+		//
+		// Decay mode is M-electron capture. In this implementation it is added to L-shell case
+		//
+		if (modeFirstRecord[5])
+		  {modeFirstRecord[5] = false; modeTotalBR[5] = b;}
+		else {
+		  G4MshellECDecayChannel *aMECChannel = new
+		    G4MshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
+					    b, c*MeV, a*MeV);
+		  aMECChannel->SetICM(applyICM);
+		  aMECChannel->SetARM(applyARM);
+		  aMECChannel->SetHLThreshold(halflifethreshold);
+		  theDecayTable->Insert(aMECChannel);
+		  modeSumBR[5] += b;
+		}
+		break;
+	      case Alpha:
+		//
+		//
+		// Decay mode is alpha.
+		//
+		if (modeFirstRecord[6])
+		  {modeFirstRecord[6] = false; modeTotalBR[6] = b;}
+		else {
+		  G4AlphaDecayChannel *anAlphaChannel = new
+		    G4AlphaDecayChannel (GetVerboseLevel(), &theParentNucleus,
+					 b, c*MeV, a*MeV);
+		  theDecayTable->Insert(anAlphaChannel);
+		  //	      delete anAlphaChannel;
+		  modeSumBR[6] += b;
+		}
+		break;
+	      case ERROR:
+	      default:
+		G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
+			    FatalException, "Error in decay mode selection");
+		
+	      }
+	    }
+	  }
+	}
+      }
+    }	
+    //
+    //
+    // Go through the decay table and make sure that the branching ratios are
+    // correctly normalised.
+    //
+    G4VDecayChannel       *theChannel             = 0;
+    G4NuclearDecayChannel *theNuclearDecayChannel = 0;
+    G4String mode                     = "";
+    G4int j                           = 0;
+    G4double theBR                    = 0.0;
+    for (i=0; i<theDecayTable->entries(); i++) {
+      theChannel             = theDecayTable->GetDecayChannel(i);
+      theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
+      theDecayMode           = theNuclearDecayChannel->GetDecayMode();
+      j          = 0;
+      if (theDecayMode != IT) {
+	theBR = theChannel->GetBR();
+	theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
+      }
+    } 
+  }		
+  DecaySchemeFile.close();		
+  if (!found && E > 0.) {
+    // cases where IT cascade for exited isotopes without entry in RDM database
+    // Decay mode is isomeric transition.
+    //
+    G4ITDecayChannel *anITChannel = new G4ITDecayChannel
+      (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
+    anITChannel->SetICM(applyICM);
+    anITChannel->SetARM(applyARM);
+    anITChannel->SetHLThreshold(halflifethreshold);
+    theDecayTable->Insert(anITChannel);
+  } 
+  if (!theDecayTable) {
     //
     // There is no radioactive decay data for this nucleus.  Return a null
     // decay table.
     //
-    {
-      G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
-      theDecayTable = 0;
-      return theDecayTable;
-    }
-  //
-  //
-  // Initialise variables used for reading in radioactive decay data.
-  //
-  G4int    nMode = 7;
-  G4bool   modeFirstRecord[7];
-  G4double modeTotalBR[7];
-  G4double modeSumBR[7];
-  G4int i;
-  for (i=0; i<nMode; i++)
-    {
-      modeFirstRecord[i] = true;
-      modeSumBR[i]       = 0.0;
-    }
-
-  G4bool complete(false);
-  G4bool found(false);
-  char inputChars[80]={' '};
-  G4String inputLine;
-  G4String recordType("");
-  G4RadioactiveDecayMode theDecayMode;
-  G4double a(0.0);
-  G4double b(0.0);
-  G4double c(0.0);
-  G4double n(1.0);
-  G4double e0;
-  //
-  //
-  // Go through each record in the data file until you identify the decay
-  // data relating to the nuclide of concern.
-  //
-//  while (!complete && -DecaySchemeFile.getline(inputChars, 80).eof() != EOF)
-  while (!complete && !DecaySchemeFile.getline(inputChars, 80).eof())    {
-      inputLine = inputChars;
-      //    G4String::stripType stripend(1);
-      //    inputLine = inputLine.strip(stripend);
-      inputLine = inputLine.strip(1);
-      if (inputChars[0] != '#' && inputLine.length() != 0)
-	{
-	  std::istringstream tmpStream(inputLine);
-	  if (inputChars[0] == 'P')
-	    //
-	    //
-	    // Nucleus is a parent type.  Check the excitation level to see if it matches
-	    // that of theParentNucleus
-	    //
-	    {
-	      tmpStream >>recordType >>a >>b;
-	      if (found) {complete = true;}
-	      else {found = (std::abs(a*keV - E)<levelTolerance);}
-	    }
-	  else if (found)
-	    {
-	      //
-	      //
-	      // The right part of the radioactive decay data file has been found.  Search
-	      // through it to determine the mode of decay of the subsequent records.
-	      //
-	      if (inputChars[0] == 'W') {
-#ifdef G4VERBOSE
-		if (GetVerboseLevel()>0) {
-		  // a comment line identified and print out the message
-		  //
-		  G4cout << " Warning in G4RadioactiveDecay::LoadDecayTable " << G4endl;
-		  G4cout << "   In data file " << file << G4endl;
-		  G4cout << "   " << inputLine << G4endl;
-		}
-#endif
-	      }	
-	      else 
-		{
-		  tmpStream >>theDecayMode >>a >>b >>c;
-		  a/=1000.;
-		  c/=1000.;
-	  
-		  //	cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl;
-	  
-		  switch (theDecayMode)
-		    {
-		    case IT:
-		      //
-		      //
-		      // Decay mode is isomeric transition.
-		      //
-		      {
-			G4ITDecayChannel *anITChannel = new G4ITDecayChannel
-			  (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, b);
-			theDecayTable->Insert(anITChannel);
-			break;
-		      }
-		    case BetaMinus:
-		      //
-		      //
-		      // Decay mode is beta-.
-		      //
-		      if (modeFirstRecord[1])
-			{modeFirstRecord[1] = false; modeTotalBR[1] = b;}
-		      else
-			{
-			  if (c > 0.) {
-			    // to work out the Fermi function normalization factor first
-			    G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, (Z+1));
-			    e0 = c*MeV/0.511;
-			    n = aBetaFermiFunction->GetFFN(e0);
-			    
-			    // now to work out the histogram and initialise the random generator
-			    G4int npti = 100;				
-			    G4double* pdf = new G4double[npti];
-			    G4int ptn;
-			    G4double g,e,ee,f;
-			    ee = e0+1.;
-			    for (ptn=0; ptn<npti; ptn++) {
-			      // e =e0*(ptn+1.)/102.;
-			      // bug fix (#662) (flei, 22/09/2004)
-			      e =e0*(ptn+0.5)/100.;
-			      g = e+1.;
-			      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
-			      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
-			    }		  
-			    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
-			    G4BetaMinusDecayChannel *aBetaMinusChannel = new
-			      G4BetaMinusDecayChannel (GetVerboseLevel(), &theParentNucleus,
-						       b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
-			    theDecayTable->Insert(aBetaMinusChannel);
-			    modeSumBR[1] += b;
-			    delete[] pdf;
-			    delete aBetaFermiFunction;
-			  }
-			}
-		      break;
-		    case BetaPlus:
-		      //
-		      //
-		      // Decay mode is beta+.
-		      //
-		      if (modeFirstRecord[2])
-			{modeFirstRecord[2] = false; modeTotalBR[2] = b;}
-		      else
-			{
-			  //  e0 = c*MeV/0.511;
-			  // bug fix (#662) (flei, 22/09/2004)
-			  // need to test e0 as there are some data files (e.g. z67.a162) which have entries for beta+
-			  // with Q < 2Me
-			  //
-			  e0 = c*MeV/0.511 -2.;
-			  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;				
-			    G4double* pdf = new G4double[npti];
-			    G4int ptn;
-			    G4double g,e,ee,f;
-			    ee = e0+1.;
-			    for (ptn=0; ptn<npti; ptn++) {
-			      // e =e0*(ptn+1.)/102.;
-			      // bug fix (#662) (flei, 22/09/2004)
-			      e =e0*(ptn+0.5)/100.;
-			      g = e+1.;
-			      f = std::sqrt(g*g-1)*(ee-g)*(ee-g)*g;
-			      pdf[ptn] = f*aBetaFermiFunction->GetFF(e);
-			    }
-			    RandGeneral* aRandomEnergy = new RandGeneral( pdf, npti);  
-			    G4BetaPlusDecayChannel *aBetaPlusChannel = new 
-			      G4BetaPlusDecayChannel (GetVerboseLevel(), &theParentNucleus,
-						      b, c*MeV, a*MeV, n, FBeta, aRandomEnergy);
-			    theDecayTable->Insert(aBetaPlusChannel);
-			    modeSumBR[2] += b;
-			    
-			    delete[] pdf;
-			    delete aBetaFermiFunction;	      
-			  }
-			}
-		      break;
-		    case KshellEC:
-		      //
-		      //
-		      // Decay mode is K-electron capture.
-		      //
-		      if (modeFirstRecord[3])
-			{modeFirstRecord[3] = false; modeTotalBR[3] = b;}
-		      else
-			{
-			  G4KshellECDecayChannel *aKECChannel = new
-			    G4KshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
-						    b, c*MeV, a*MeV);
-			  theDecayTable->Insert(aKECChannel);
-			  //delete aKECChannel;
-			  modeSumBR[3] += b;
-			}
-		      break;
-		    case LshellEC:
-		      //
-		      //
-		      // Decay mode is L-electron capture.
-		      //
-		      if (modeFirstRecord[4])
-			{modeFirstRecord[4] = false; modeTotalBR[4] = b;}
-		      else
-			{
-			  G4LshellECDecayChannel *aLECChannel = new
-			    G4LshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
-						    b, c*MeV, a*MeV);
-			  theDecayTable->Insert(aLECChannel);
-			  //delete aLECChannel;
-			  modeSumBR[4] += b;
-			}
-		      break;
-		    case MshellEC:
-		      //
-		      //
-		      // Decay mode is M-electron capture. In this implementation it is added to L-shell case
-		      //
-		      if (modeFirstRecord[5])
-			{modeFirstRecord[5] = false; modeTotalBR[5] = b;}
-		      else
-			{
-			  G4MshellECDecayChannel *aMECChannel = new
-			    G4MshellECDecayChannel (GetVerboseLevel(), &theParentNucleus,
-						    b, c*MeV, a*MeV);
-			  theDecayTable->Insert(aMECChannel);
-			  modeSumBR[5] += b;
-			}
-		      break;
-		    case Alpha:
-		      //
-		      //
-		      // Decay mode is alpha.
-		      //
-		      if (modeFirstRecord[6])
-			{modeFirstRecord[6] = false; modeTotalBR[6] = b;}
-		      else
-			{
-			  G4AlphaDecayChannel *anAlphaChannel = new
-			    G4AlphaDecayChannel (GetVerboseLevel(), &theParentNucleus,
-						 b, c*MeV, a*MeV);
-			  theDecayTable->Insert(anAlphaChannel);
-			  //	      delete anAlphaChannel;
-			  modeSumBR[6] += b;
-			}
-		      break;
-		    case ERROR:
-		    default:
-                      G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601",
-                                  FatalException, "Error in decay mode selection");
-		      
-		    }
-		}
-	    }
-	}
-
-    }  
-  DecaySchemeFile.close();
-  if (!found && E > 0.) {
-    // cases where IT cascade follows immediately after a decay
-    
-    // Decay mode is isomeric transition.
-    //
-    G4ITDecayChannel *anITChannel = new G4ITDecayChannel
-      (GetVerboseLevel(), (const G4Ions*) &theParentNucleus, 1);
-    theDecayTable->Insert(anITChannel);
-  }
-  //
-  //
-  // Go through the decay table and make sure that the branching ratios are
-  // correctly normalised.
-  //
-  G4VDecayChannel       *theChannel             = 0;
-  G4NuclearDecayChannel *theNuclearDecayChannel = 0;
-  G4String mode                     = "";
-  G4int j                           = 0;
-  G4double theBR                    = 0.0;
-  for (i=0; i<theDecayTable->entries(); i++)
-    {
-      theChannel             = theDecayTable->GetDecayChannel(i);
-      theNuclearDecayChannel = static_cast<G4NuclearDecayChannel *>(theChannel);
-      theDecayMode           = theNuclearDecayChannel->GetDecayMode();
-      j          = 0;
-      if (theDecayMode != IT)
-	{
-	  theBR = theChannel->GetBR();
-	  theChannel->SetBR(theBR*modeTotalBR[theDecayMode]/modeSumBR[theDecayMode]);
-	}
-    }  
-
+    G4cerr <<"G4RadoactiveDecay::LoadDecayTable() : cannot find ion radioactive decay file " <<G4endl;
+    theDecayTable = 0;
+    return theDecayTable;
+  }	
   if (theDecayTable && GetVerboseLevel()>1)
     {
@@ -1086,6 +1091,6 @@
 	  if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) {
 	    
-	    // Level hafe life is in ns and I want to set the gate as 1 micros
-	    if ( theDecayMode == 0 && level->HalfLife() >= 1000.){    
+	    // 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
 	      theDecayTable->Insert(theChannel);
@@ -1424,4 +1429,14 @@
 	G4DecayProducts *products = DoDecay(*theParticleDef);
 	//
+	// check if the product is the same as input and kill the track if necessary to prevent infinite loop
+	// (11/05/10, F.Lei)
+	//
+	if ( products->entries() == 1) {
+	  fParticleChangeForRadDecay.SetNumberOfSecondaries(0);
+	  fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ;
+	  fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0);
+	  ClearNumberOfInteractionLengthLeft();
+	  return &fParticleChangeForRadDecay;
+	}	
 	//
 	// Get parent particle information and boost the decay products to the
@@ -1440,6 +1455,8 @@
 	    // time lapsed between the particle come to rest and the actual decay. This time 
 	    // is simply sampled with the mean-life of the particle.
-	    //
-	    finalGlobalTime += -std::log( G4UniformRand()) * theParticleDef->GetPDGLifeTime() ;
+	    // but we need to protect the case PDGTime < 0. (F.Lei 11/05/10)
+	    G4double temptime =  -std::log( G4UniformRand()) * theParticleDef->GetPDGLifeTime();
+	    if (temptime <0.) temptime =0.; 
+	    finalGlobalTime += temptime ;
 	    energyDeposit += theParticle->GetKineticEnergy();
 	  }
Index: trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecaymessenger.cc
===================================================================
--- trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecaymessenger.cc	(revision 1201)
+++ trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecaymessenger.cc	(revision 1315)
@@ -116,5 +116,27 @@
   brbiasCmd->SetParameterName("BRBias",true);
   brbiasCmd->SetDefaultValue(true);
-
+  //
+  // Command contols whether ICM will be applied or not
+  //
+  icmCmd = new G4UIcmdWithABool ("/grdm/applyICM",this);
+  icmCmd->SetGuidance("True: ICM is applied; false: no");
+  icmCmd->SetParameterName("applyICM",true);
+  icmCmd->SetDefaultValue(true);
+  //
+  // Command contols whether ARM will be applied or not
+  //
+  armCmd = new G4UIcmdWithABool ("/grdm/applyARM",this);
+  armCmd->SetGuidance("True: ARM is applied; false: no");
+  armCmd->SetParameterName("applyARM",true);
+  armCmd->SetDefaultValue(true);
+  //
+  // Command to set the h-l thresold for isomer production
+  //
+  hlthCmd = new G4UIcmdWithADoubleAndUnit("/grdm/hlThreshold",this);
+  hlthCmd->SetGuidance("Set the h-l threshold for isomer production");
+  hlthCmd->SetParameterName("hlThreshold",false);
+  hlthCmd->SetRange("hlThreshold>0.");
+  hlthCmd->SetUnitCategory("Time");
+  hlthCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
   //
   // Command to define the incident particle source time profile.
@@ -174,4 +196,7 @@
   delete allvolumesCmd;
   delete deallvolumesCmd;
+  delete icmCmd;
+  delete armCmd;
+  delete hlthCmd;
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -215,4 +240,10 @@
   else if (command==verboseCmd) {theRadioactiveDecayContainer->
 				       SetVerboseLevel(verboseCmd->GetNewIntValue(newValues));}
+  else if (command==icmCmd ) {theRadioactiveDecayContainer->
+      SetICM(icmCmd->GetNewBoolValue(newValues));}
+  else if (command==armCmd ) {theRadioactiveDecayContainer->
+      SetARM(armCmd->GetNewBoolValue(newValues));}
+  else if (command==hlthCmd ) {theRadioactiveDecayContainer->
+      SetHLThreshold(hlthCmd->GetNewDoubleValue(newValues));}
 }
 
