Index: trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationRudd.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationRudd.cc	(revision 819)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4FinalStateIonisationRudd.cc	(revision 961)
@@ -24,53 +24,13 @@
 // ********************************************************************
 //
-//
-// $Id: G4FinalStateIonisationRudd.cc,v 1.5 2007/11/26 17:27:09 pia Exp $
-// GEANT4 tag $Name:  $
-// 
-// Contact Author: Sebastien Incerti (incerti@cenbg.in2p3.fr)
-//                 Maria Grazia Pia  (Maria.Grazia.Pia@cern.ch)
-//
-///
-// Reference: TNS Geant4-DNA paper
-// Reference for implementation model: NIM. 155, pp. 145-156, 1978
-//
-// History:
-// -----------
-// Date         Name              Modification
-// 28 Apr 2007  M.G. Pia          Created in compliance with design described in TNS paper
-//    Nov 2007  S. Incerti        Implementation
-// 26 Nov 2007  MGP               Cleaned up std::
-//
-// -------------------------------------------------------------------
-
-// Class description:
-// Reference: TNS Geant4-DNA paper
-// S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
-// design foundation and implementation of the first set of models,
-// IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
-// Further documentation available from http://www.ge.infn.it/geant4/dna
-
-// -------------------------------------------------------------------
-
+// $Id: G4FinalStateIonisationRudd.cc,v 1.8 2008/08/20 14:51:48 sincerti Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 
 #include "G4FinalStateIonisationRudd.hh"
-#include "G4Track.hh"
-#include "G4Step.hh"
-#include "G4DynamicParticle.hh"
-#include "Randomize.hh"
-
-#include "G4ParticleTypes.hh"
-#include "G4ParticleDefinition.hh"
-#include "G4Electron.hh"
-#include "G4Proton.hh"
-#include "G4SystemOfUnits.hh"
-#include "G4ParticleMomentum.hh"
-#include "G4DNAGenericIonsManager.hh"
-
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4FinalStateIonisationRudd::G4FinalStateIonisationRudd()
 {
-  name = "IonisationBorn";
-  // Default energy limits (defined for protection against anomalous behaviour only)
   lowEnergyLimitDefault = 100 * eV;
   highEnergyLimitDefault = 100 * MeV;
@@ -112,13 +72,13 @@
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4FinalStateIonisationRudd::~G4FinalStateIonisationRudd()
-{ }
-
-
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 const G4FinalStateProduct& G4FinalStateIonisationRudd::GenerateFinalState(const G4Track& track, const G4Step& /* step */)
 {
-  // Clear previous secondaries, energy deposit and particle kill status
   product.Clear();
 
@@ -132,30 +92,22 @@
   const G4String& particleName = particle->GetDefinition()->GetParticleName();
 
-  // Retrieve energy limits for the current particle type
-
   std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
   pos1 = lowEnergyLimit.find(particleName);
 
-  // Lower limit
   if (pos1 != lowEnergyLimit.end())
-    {
-      lowLim = pos1->second;
-    }
-
-  // Upper limit
+  {
+    lowLim = pos1->second;
+  }
+
   std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
   pos2 = highEnergyLimit.find(particleName);
 
   if (pos2 != highEnergyLimit.end())
-    {
-      highLim = pos2->second;
-    }
-
-  // Verify that the current track is within the energy limits of validity of the cross section model
+  {
+    highLim = pos2->second;
+  }
 
   if (k >= lowLim && k <= highLim)
-    {
-      // Kinetic energy of primary particle
-
+  {
       G4ParticleDefinition* definition = particle->GetDefinition();
       G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
@@ -186,5 +138,4 @@
       G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
 
-      // Primary Particle Direction
       G4double finalPx = totalMomentum*primaryDirection.x() - deltaTotalMomentum*deltaDirection.x();
       G4double finalPy = totalMomentum*primaryDirection.y() - deltaTotalMomentum*deltaDirection.y();
@@ -200,12 +151,15 @@
       G4DynamicParticle* aElectron = new G4DynamicParticle(G4Electron::Electron(),deltaDirection,secondaryKinetic);
       product.AddSecondary(aElectron);
-    }
-
-  if (k < lowLim) {product.KillPrimaryParticle();product.AddEnergyDeposit(k);}
+  }
+
+  if (k < lowLim) 
+  {  
+    product.KillPrimaryParticle();
+  }
  
   return product;
 }
 
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4double G4FinalStateIonisationRudd::RandomizeEjectedElectronEnergy(G4ParticleDefinition* particleDefinition, 
@@ -220,24 +174,27 @@
   if (particleDefinition == G4Proton::ProtonDefinition() 
       || particleDefinition == instance->GetIon("hydrogen"))
-	  
-    { 
+  { 
       maximumKineticEnergyTransfer= 4.* (electron_mass_c2 / proton_mass_c2) * k;
-    }
+  }
 
   if (particleDefinition == instance->GetIon("helium") 
       || particleDefinition == instance->GetIon("alpha+")
       || particleDefinition == instance->GetIon("alpha++"))
-    { 
+  { 
       maximumKineticEnergyTransfer= 4.* (0.511 / 3728) * k;
-    }
+  }
 
   G4double crossSectionMaximum = 0.;
+  
   for(G4double value=waterStructure.IonisationEnergy(shell); value<=4.*waterStructure.IonisationEnergy(shell) ; value+=0.1*eV)
-    {
+  {
       G4double differentialCrossSection = DifferentialCrossSection(particleDefinition, k, value, shell);
       if(differentialCrossSection >= crossSectionMaximum) crossSectionMaximum = differentialCrossSection;
-    }
+  }
+  
   G4double secElecKinetic = 0.;
-  do{
+  
+  do
+  {
     secElecKinetic = G4UniformRand() * maximumKineticEnergyTransfer;
   } while(G4UniformRand()*crossSectionMaximum > DifferentialCrossSection(particleDefinition, 
@@ -249,10 +206,12 @@
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
 
 void G4FinalStateIonisationRudd::RandomizeEjectedElectronDirection(G4ParticleDefinition* particleDefinition, 
 								   G4double k, 
 								   G4double secKinetic, 
-								   G4double cosTheta, 
-								   G4double phi )
+								   G4double & cosTheta, 
+								   G4double & phi )
 {
   G4DNAGenericIonsManager *instance;
@@ -263,18 +222,20 @@
   if (particleDefinition == G4Proton::ProtonDefinition() 
       || particleDefinition == instance->GetIon("hydrogen")) 
-    { 
+  { 
       maxSecKinetic = 4.* (electron_mass_c2 / proton_mass_c2) * k;
-    }
+  }
   
   if (particleDefinition == instance->GetIon("helium") 
       || particleDefinition == instance->GetIon("alpha+")
       || particleDefinition == instance->GetIon("alpha++"))
-    { 
+  { 
       maxSecKinetic = 4.* (0.511 / 3728) * k;
-    }
+  }
   
   phi = twopi * G4UniformRand();
   cosTheta = std::sqrt(secKinetic / maxSecKinetic);
 }
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 
@@ -314,5 +275,5 @@
 
   if (j == 4) 
-    {
+  {
       //Data For Liquid Water K SHELL from Dingfelder (Protons in Water)
       A1 = 1.25; 
@@ -326,7 +287,7 @@
       D2 = 0.00; 
       alphaConst = 0.66;
-    }
+  }
   else 
-    {
+  {
       //Data For Liquid Water from Dingfelder (Protons in Water)
       A1 = 1.02; 
@@ -340,11 +301,8 @@
       D2 = 0.04; 
       alphaConst = 0.64;
-    }
+  }
   
   const G4double n = 2.;
   const G4double Gj[5] = {0.99, 1.11, 1.11, 0.52, 1.};
-
-  //const G4double I[5]={12.61*eV, 14.73*eV, 18.55*eV, 32.2*eV, 539.7*eV}; // for water Vapor
-  //const G4double energyConstant[]={10.79*eV, 13.39*eV, 16.05*eV, 32.30*eV, 539.*eV};
 
   G4DNAGenericIonsManager* instance;
@@ -359,14 +317,14 @@
   if (particleDefinition == G4Proton::ProtonDefinition() 
       || particleDefinition == instance->GetIon("hydrogen")) 
-    {
+  {
       tau = (electron_mass_c2/proton_mass_c2) * k ;
-    }
+  }
    
   if ( particleDefinition == instance->GetIon("helium") 
        || particleDefinition == instance->GetIon("alpha+")
        || particleDefinition == instance->GetIon("alpha++")) 
-    {
+  {
       tau = (0.511/3728.) * k ;
-    }
+  }
  
   G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
@@ -390,12 +348,10 @@
 	  || particleDefinition == instance->GetIon("hydrogen")
 	  ) 
-    {
+  {
       return(sigma);
-    }
-
-  // ------------
-  
+  }
+
   if (particleDefinition == instance->GetIon("alpha++") ) 
-    {
+  {
       slaterEffectiveCharge[0]=0.;
       slaterEffectiveCharge[1]=0.;
@@ -404,8 +360,8 @@
       sCoefficient[1]=0.;
       sCoefficient[2]=0.;
-    }
+  }
 
   if (particleDefinition == instance->GetIon("alpha+") ) 
-    {
+  {
       slaterEffectiveCharge[0]=2.0;
       slaterEffectiveCharge[1]=1.15;
@@ -414,8 +370,8 @@
       sCoefficient[1]=0.15;
       sCoefficient[2]=0.15;
-    }
+  }
 
   if (particleDefinition == instance->GetIon("helium") ) 
-    {
+  {
       slaterEffectiveCharge[0]=1.7;
       slaterEffectiveCharge[1]=1.15;
@@ -424,5 +380,5 @@
       sCoefficient[1]=0.25;
       sCoefficient[2]=0.25;
-    }
+  }
   
   if (    particleDefinition == instance->GetIon("helium") 
@@ -430,5 +386,5 @@
 	  || particleDefinition == instance->GetIon("alpha++")
 	  ) 
-    {
+  {
       sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
     
@@ -440,8 +396,10 @@
 	   
       return zEff * zEff * sigma ;
-    }  
+  }  
   
   return 0;
 }
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4double G4FinalStateIonisationRudd::S_1s(G4double t, 
@@ -459,5 +417,5 @@
 }
 
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4double G4FinalStateIonisationRudd::S_2s(G4double t,
@@ -476,5 +434,5 @@
 }
 
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4double G4FinalStateIonisationRudd::S_2p(G4double t, 
@@ -492,5 +450,5 @@
 }
 
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4double G4FinalStateIonisationRudd::R(G4double t,
@@ -500,5 +458,4 @@
 {
   // tElectron = m_electron / m_alpha * t
-  // Hardcoded in Riccardo's implementation; to be corrected
   // Dingfelder, in Chattanooga 2005 proceedings, p 4
 
@@ -509,5 +466,5 @@
 }
 
-
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4double G4FinalStateIonisationRudd::CorrectionFactor(G4ParticleDefinition* particleDefinition, G4double k) 
@@ -517,16 +474,16 @@
 
   if (particleDefinition == G4Proton::Proton()) 
-    {
+  {
       return(1.);
-    }
+  }
   else 
     if (particleDefinition == instance->GetIon("hydrogen")) 
-      { 
+    { 
 	G4double value = (std::log(k/eV)-4.2)/0.5;
 	return((0.8/(1+std::exp(value))) + 0.9);
-      }
+    }
     else 
-      {    
+    {    
 	return(1.);
-      }
-}
+    }
+}
