Index: trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationModel.cc	(revision 1337)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4DNARuddIonisationModel.cc	(revision 1340)
@@ -24,6 +24,6 @@
 // ********************************************************************
 //
-// $Id: G4DNARuddIonisationModel.cc,v 1.17 2010/04/07 20:08:31 sincerti Exp $
-// GEANT4 tag $Name: geant4-09-04-beta-01 $
+// $Id: G4DNARuddIonisationModel.cc,v 1.21 2010/11/04 14:52:17 sincerti Exp $
+// GEANT4 tag $Name: emlowen-V09-03-54 $
 //
 
@@ -159,5 +159,5 @@
 
     lowEnergyLimit[alphaPlusPlus] = lowEnergyLimitForZ2;
-    highEnergyLimit[alphaPlusPlus] = 10. * MeV;
+    highEnergyLimit[alphaPlusPlus] = 400. * MeV;
 
     // Cross section
@@ -181,5 +181,5 @@
 
     lowEnergyLimit[alphaPlus] = lowEnergyLimitForZ2;
-    highEnergyLimit[alphaPlus] = 10. * MeV;
+    highEnergyLimit[alphaPlus] = 400. * MeV;
 
     // Cross section
@@ -202,5 +202,5 @@
 
     lowEnergyLimit[helium] = lowEnergyLimitForZ2;
-    highEnergyLimit[helium] = 10. * MeV;
+    highEnergyLimit[helium] = 400. * MeV;
 
     // Cross section
@@ -363,50 +363,4 @@
 	      sigma = table->FindValue(k);
 
-	      // BEGIN ELECTRON CORRECTION
-	      // add ONE or TWO electron-water excitation for alpha+ and helium
-   
-	      if ( particleDefinition == instance->GetIon("alpha+") 
-		   ||
-		   particleDefinition == instance->GetIon("helium")
-		   ) 
-	      {
-      
-		  G4DNACrossSectionDataSet* electronDataset = new G4DNACrossSectionDataSet 
-		    (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
-       
-		  electronDataset->LoadData("dna/sigma_ionisation_e_born");
-
-		  G4double kElectron = k * 0.511/3728;
-
-		  if ( particleDefinition == instance->GetIon("alpha+") ) 
-		  {
-		      G4double tmp1 = table->FindValue(k) + electronDataset->FindValue(kElectron);
-		      delete electronDataset;
-		      if (verboseLevel > 3)
-                      {
-                        G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
-                        G4cout << " - Cross section per water molecule (cm^2)=" << tmp1/cm/cm << G4endl;
-                        G4cout << " - Cross section per water molecule (cm^-1)=" << 
-			tmp1*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
-                      } 
-		      return tmp1*material->GetAtomicNumDensityVector()[1];
-		  }
-
-		  if ( particleDefinition == instance->GetIon("helium") ) 
-		  {
-		      G4double tmp2 = table->FindValue(k) +  2. * electronDataset->FindValue(kElectron);
-		      delete electronDataset;
-		      if (verboseLevel > 3)
-                      {
-                        G4cout << "---> Kinetic energy(eV)=" << k/eV << G4endl;
-                        G4cout << " - Cross section per water molecule (cm^2)=" << tmp2/cm/cm << G4endl;
-                        G4cout << " - Cross section per water molecule (cm^-1)=" << tmp2*
-			material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
-                      } 
-		      return tmp2*material->GetAtomicNumDensityVector()[1];
-		  }
-              }      
-
-	      // END ELECTRON CORRECTION
          }
       }
@@ -489,8 +443,10 @@
       G4ParticleDefinition* definition = particle->GetDefinition();
       G4ParticleMomentum primaryDirection = particle->GetMomentumDirection();
+      /*
       G4double particleMass = definition->GetPDGMass();
       G4double totalEnergy = k + particleMass;
       G4double pSquare = k*(totalEnergy+particleMass);
       G4double totalMomentum = std::sqrt(pSquare);
+      */
 
       G4int ionizationShell = RandomSelect(k,particleName);
@@ -511,4 +467,6 @@
       deltaDirection.rotateUz(primaryDirection);
 
+      // Ignored for ions on electrons
+      /*
       G4double deltaTotalMomentum = std::sqrt(secondaryKinetic*(secondaryKinetic + 2.*electron_mass_c2 ));
 
@@ -525,4 +483,7 @@
 
       fParticleChangeForGamma->ProposeMomentumDirection(direction.unit()) ;
+      */
+      fParticleChangeForGamma->ProposeMomentumDirection(primaryDirection);
+
       fParticleChangeForGamma->SetProposedKineticEnergy(k-bindingEnergy-secondaryKinetic);
       fParticleChangeForGamma->ProposeLocalEnergyDeposit(bindingEnergy);
@@ -575,4 +536,5 @@
   }
   
+ 
   G4double secElecKinetic = 0.;
   
@@ -616,5 +578,12 @@
   
   phi = twopi * G4UniformRand();
-  cosTheta = std::sqrt(secKinetic / maxSecKinetic);
+  
+  //cosTheta = std::sqrt(secKinetic / maxSecKinetic);
+
+  // Restriction below 100 eV from Emfietzoglou (2000)
+
+  if (secKinetic>100*eV) cosTheta = std::sqrt(secKinetic / maxSecKinetic);
+  else cosTheta = (2.*G4UniformRand())-1.;
+  
 }
 
@@ -656,5 +625,7 @@
   G4double alphaConst ;
 
-  const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
+//  const G4double Bj[5] = {12.61*eV, 14.73*eV, 18.55*eV, 32.20*eV, 539.7*eV};
+// The following values are provided by M. dingfelder (priv. comm)
+  const G4double Bj[5] = {12.60*eV, 14.70*eV, 18.40*eV, 32.20*eV, 540*eV};
 
   if (j == 4) 
@@ -681,5 +652,7 @@
       E1 = 0.38; 
       A2 = 1.07; 
-      B2 = 14.6;
+      // Value provided by M. Dingfelder (priv. comm)
+      B2 = 11.6;
+      //
       C2 = 0.60; 
       D2 = 0.04; 
@@ -697,4 +670,7 @@
 
   G4double w = wBig / Bj[ionizationLevelIndex];
+  // Note that the following (j==4) cases are provided by   M. Dingfelder (priv. comm)
+  if (j==4) w = wBig / waterStructure.IonisationEnergy(ionizationLevelIndex);
+
   G4double Ry = 13.6*eV;
 
@@ -713,9 +689,14 @@
       tau = (0.511/3728.) * k ;
   }
- 
+
   G4double S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/Bj[ionizationLevelIndex]),2);
+  if (j==4) S = 4.*pi * Bohr_radius*Bohr_radius * n * std::pow((Ry/waterStructure.IonisationEnergy(ionizationLevelIndex)),2);
+
   G4double v2 = tau / Bj[ionizationLevelIndex];
+  if (j==4) v2 = tau / waterStructure.IonisationEnergy(ionizationLevelIndex);
+
   G4double v = std::sqrt(v2);
   G4double wc = 4.*v2 - 2.*v - (Ry/(4.*Bj[ionizationLevelIndex]));
+  if (j==4) wc = 4.*v2 - 2.*v - (Ry/(4.*waterStructure.IonisationEnergy(ionizationLevelIndex)));
 
   G4double L1 = (C1* std::pow(v,(D1))) / (1.+ E1*std::pow(v, (D1+4.)));
@@ -731,7 +712,12 @@
     * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
 
+  if (j==4) sigma = CorrectionFactor(particleDefinition, k) 
+    * Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) 
+    * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
+
   if ( (particleDefinition == instance->GetIon("hydrogen")) && (ionizationLevelIndex==4)) 
 
-    sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) 
+//    sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) 
+    sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) 
     * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
 
@@ -756,6 +742,8 @@
   {
       slaterEffectiveCharge[0]=2.0;
-      slaterEffectiveCharge[1]=1.15;
-      slaterEffectiveCharge[2]=1.15;
+      // The following values are provided by M. Dingfelder (priv. comm)    
+      slaterEffectiveCharge[1]=2.0;
+      slaterEffectiveCharge[2]=2.0;
+      //
       sCoefficient[0]=0.7;
       sCoefficient[1]=0.15;
@@ -780,4 +768,7 @@
       sigma = Gj[j] * (S/Bj[ionizationLevelIndex]) * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
     
+      if (j==4) sigma = Gj[j] * (S/waterStructure.IonisationEnergy(ionizationLevelIndex)) 
+                              * ( (F1+w*F2) / ( std::pow((1.+w),3) * ( 1.+std::exp(alphaConst*(w-wc)/v))) );
+
       G4double zEff = particleDefinition->GetPDGCharge() / eplus + particleDefinition->GetLeptonNumber();
   
@@ -852,5 +843,7 @@
 
   G4double tElectron = 0.511/3728. * t;
-  G4double value = 2. * tElectron * slaterEffectiveChg / (energyTransferred * shellNumber);
+  // The following values are provided by M. Dingfelder (priv. comm)    
+  G4double H = 2.*13.60569172 * eV;
+  G4double value = std::sqrt ( 2. * tElectron / H ) / ( energyTransferred / H ) *  (slaterEffectiveChg/shellNumber);
   
   return value;
@@ -872,5 +865,6 @@
     { 
 	G4double value = (std::log10(k/eV)-4.2)/0.5;
-	return((0.8/(1+std::exp(value))) + 0.9);
+	// The following values are provided by M. Dingfelder (priv. comm)    
+        return((0.6/(1+std::exp(value))) + 0.9);
     }
     else 
@@ -891,21 +885,5 @@
   G4DNAGenericIonsManager *instance;
   instance = G4DNAGenericIonsManager::Instance();
-  G4double kElectron(0);
-
-  G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
  
-  if ( particle == instance->GetIon("alpha+")->GetParticleName()
-       ||
-       particle == instance->GetIon("helium")->GetParticleName()
-       ) 
-  {     
-      electronDataset->LoadData("dna/sigma_ionisation_e_born");
-
-      kElectron = k * 0.511/3728;
-       
-  }      
-  
-  // END PART 1/2 OF ELECTRON CORRECTION
-  
   G4int level = 0;
 
@@ -931,16 +909,4 @@
 	      i--;
 	      valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
-
-	      // BEGIN PART 2/2 OF ELECTRON CORRECTION
-	      // Use only electron partial cross sections
-	      
-	      if (particle == instance->GetIon("alpha+")->GetParticleName()) 
-		{valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronDataset->GetComponent(i)->FindValue(kElectron); }
-
-	      if (particle == instance->GetIon("helium")->GetParticleName()) 
-		{valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronDataset->GetComponent(i)->FindValue(kElectron); }
-
-	      // BEGIN PART 2/2 OF ELECTRON CORRECTION
-
 	      value += valuesBuffer[i];
 	  }
@@ -958,7 +924,4 @@
 	      {
 		  delete[] valuesBuffer;
-		  
-		  if (electronDataset) delete electronDataset;
-		  
 		  return i;
 	      }
@@ -975,6 +938,4 @@
   }
 
-  delete electronDataset;
-      
   return level;
 }
