Index: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedGammaConversionModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedGammaConversionModel.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedGammaConversionModel.cc	(revision 1197)
@@ -0,0 +1,945 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+
+#include "G4LivermorePolarizedGammaConversionModel.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+using namespace std;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4LivermorePolarizedGammaConversionModel::G4LivermorePolarizedGammaConversionModel(const G4ParticleDefinition*,
+                                             const G4String& nam)
+  :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),crossSectionHandler(0)
+{
+  lowEnergyLimit = 1.0220000 * MeV; 
+  highEnergyLimit = 100 * GeV;
+  SetLowEnergyLimit(lowEnergyLimit);
+  SetHighEnergyLimit(highEnergyLimit);
+  smallEnergy = 2.*MeV;
+
+
+  verboseLevel= 0;
+  // Verbosity scale:
+  // 0 = nothing 
+  // 1 = warning for energy non-conservation 
+  // 2 = details of energy budget
+  // 3 = calculation of cross sections, file openings, samping of atoms
+  // 4 = entering in methods
+
+  G4cout << "Livermore Polarized GammaConversion is constructed " << G4endl
+         << "Energy range: "
+         << lowEnergyLimit / keV << " keV - "
+         << highEnergyLimit / GeV << " GeV"
+         << G4endl;
+
+  crossSectionHandler = new G4CrossSectionHandler();
+  crossSectionHandler->Initialise(0,1.0220*MeV,100.*GeV,400);
+
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4LivermorePolarizedGammaConversionModel::~G4LivermorePolarizedGammaConversionModel()
+{  
+  //  if (meanFreePathTable)   delete meanFreePathTable;
+  if (crossSectionHandler) delete crossSectionHandler;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LivermorePolarizedGammaConversionModel::Initialise(const G4ParticleDefinition* particle,
+                                       const G4DataVector& cuts)
+{
+  if (verboseLevel > 3)
+    G4cout << "Calling G4LivermorePolarizedGammaConversionModel::Initialise()" << G4endl;
+
+  if (crossSectionHandler)
+  {
+    crossSectionHandler->Clear();
+    delete crossSectionHandler;
+  }
+
+
+  // Energy limits
+  
+  if (LowEnergyLimit() < lowEnergyLimit)
+    {
+      G4cout << "G4LivermorePolarizedGammaConversionModel: low energy limit increased from " << 
+	LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
+      SetLowEnergyLimit(lowEnergyLimit);
+    }
+
+  if (HighEnergyLimit() > highEnergyLimit)
+    {
+      G4cout << "G4LivermorePolarizedGammaConversionModel: high energy limit decreased from " << 
+	HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
+      SetHighEnergyLimit(highEnergyLimit);
+    }
+
+  // Reading of data files - all materials are read
+  
+  crossSectionHandler = new G4CrossSectionHandler;
+  crossSectionHandler->Clear();
+  G4String crossSectionFile = "pair/pp-cs-";
+  crossSectionHandler->LoadData(crossSectionFile);
+
+  //  meanFreePathTable = 0;
+  //meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
+
+
+  //
+  if (verboseLevel > 2) 
+    G4cout << "Loaded cross section files for Livermore Polarized GammaConversion model" << G4endl;
+
+  InitialiseElementSelectors(particle,cuts);
+
+  G4cout << "Livermore Polarized GammaConversion model is initialized " << G4endl
+         << "Energy range: "
+         << LowEnergyLimit() / keV << " keV - "
+         << HighEnergyLimit() / GeV << " GeV"
+         << G4endl;
+
+  //
+    
+  if(isInitialised) return;
+
+  if(pParticleChange)
+    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
+  else
+    fParticleChange = new G4ParticleChangeForGamma();
+    
+  isInitialised = true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4LivermorePolarizedGammaConversionModel::ComputeCrossSectionPerAtom(
+                                       const G4ParticleDefinition*,
+                                             G4double GammaEnergy,
+                                             G4double Z, G4double,
+                                             G4double, G4double)
+{
+  if (verboseLevel > 3)
+    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedGammaConversionModel" << G4endl;
+
+  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
+  return cs;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LivermorePolarizedGammaConversionModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
+							       const G4MaterialCutsCouple* couple,
+							       const G4DynamicParticle* aDynamicGamma,
+							       G4double,
+							       G4double)
+{
+
+  // Fluorescence generated according to:
+  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
+  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
+  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
+
+  if (verboseLevel > 3)
+    G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedGammaConversionModel" << G4endl;
+
+  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
+  // Within energy limit?
+
+  if(photonEnergy <= lowEnergyLimit)
+    {
+      fParticleChange->ProposeTrackStatus(fStopAndKill);
+      fParticleChange->SetProposedKineticEnergy(0.);
+      return;
+    }
+
+
+  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
+  G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
+
+  // Make sure that the polarization vector is perpendicular to the
+  // gamma direction. If not
+
+  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
+    { // only for testing now
+      gammaPolarization0 = GetRandomPolarization(gammaDirection0);
+    }
+  else
+    {
+      if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
+	{
+	  gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
+	}
+    }
+
+  // End of Protection
+
+
+  G4double epsilon ;
+  G4double epsilon0 = electron_mass_c2 / photonEnergy ;
+
+  // Do it fast if photon energy < 2. MeV
+
+  if (photonEnergy < smallEnergy )
+    {
+      epsilon = epsilon0 + (0.5 - epsilon0) * G4UniformRand();
+    }
+  else
+    {
+
+ // Select randomly one element in the current material
+
+      //     G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
+      const G4Element* element = crossSectionHandler->SelectRandomElement(couple,photonEnergy);
+
+      if (element == 0)
+        {
+          G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - element = 0" << G4endl;
+        }
+      G4IonisParamElm* ionisation = element->GetIonisation();
+      if (ionisation == 0)
+        {
+          G4cout << "G4LivermorePolarizedGammaConversionModel::PostStepDoIt - ionisation = 0" << G4endl;
+        }
+
+
+
+      // Extract Coulomb factor for this Element
+
+      G4double fZ = 8. * (ionisation->GetlogZ3());
+      if (photonEnergy > 50. * MeV) fZ += 8. * (element->GetfCoulomb());
+
+      // Limits of the screening variable
+      G4double screenFactor = 136. * epsilon0 / (element->GetIonisation()->GetZ3()) ;
+      G4double screenMax = exp ((42.24 - fZ)/8.368) - 0.952 ;
+      G4double screenMin = std::min(4.*screenFactor,screenMax) ;
+
+      // Limits of the energy sampling
+      G4double epsilon1 = 0.5 - 0.5 * sqrt(1. - screenMin / screenMax) ;
+      G4double epsilonMin = std::max(epsilon0,epsilon1);
+      G4double epsilonRange = 0.5 - epsilonMin ;
+
+      // Sample the energy rate of the created electron (or positron)
+      G4double screen;
+      G4double gReject ;
+
+      G4double f10 = ScreenFunction1(screenMin) - fZ;
+      G4double f20 = ScreenFunction2(screenMin) - fZ;
+      G4double normF1 = std::max(f10 * epsilonRange * epsilonRange,0.);
+      G4double normF2 = std::max(1.5 * f20,0.);
+
+      do {
+        if (normF1 / (normF1 + normF2) > G4UniformRand() )
+          {
+            epsilon = 0.5 - epsilonRange * pow(G4UniformRand(), 0.3333) ;
+            screen = screenFactor / (epsilon * (1. - epsilon));
+            gReject = (ScreenFunction1(screen) - fZ) / f10 ;
+          }
+        else
+          {
+            epsilon = epsilonMin + epsilonRange * G4UniformRand();
+            screen = screenFactor / (epsilon * (1 - epsilon));
+            gReject = (ScreenFunction2(screen) - fZ) / f20 ;
+
+
+	  }
+      } while ( gReject < G4UniformRand() );
+
+    }   //  End of epsilon sampling
+  
+  // Fix charges randomly
+  
+  G4double electronTotEnergy;
+  G4double positronTotEnergy;
+
+
+  if (CLHEP::RandBit::shootBit())
+    {
+      electronTotEnergy = (1. - epsilon) * photonEnergy;
+      positronTotEnergy = epsilon * photonEnergy;
+    }
+  else
+    {
+      positronTotEnergy = (1. - epsilon) * photonEnergy;
+      electronTotEnergy = epsilon * photonEnergy;
+    }
+
+  // Scattered electron (positron) angles. ( Z - axis along the parent photon)
+  // Universal distribution suggested by L. Urban (Geant3 manual (1993) Phys211),
+  // derived from Tsai distribution (Rev. Mod. Phys. 49, 421 (1977)
+
+  G4double u;
+  const G4double a1 = 0.625;
+  G4double a2 = 3. * a1;
+
+  if (0.25 > G4UniformRand())
+    {
+      u = - log(G4UniformRand() * G4UniformRand()) / a1 ;
+    }
+  else
+    {
+      u = - log(G4UniformRand() * G4UniformRand()) / a2 ;
+    }
+
+  G4double Ene = electronTotEnergy/electron_mass_c2; // Normalized energy
+
+  G4double cosTheta = 0.;
+  G4double sinTheta = 0.;
+
+  SetTheta(&cosTheta,&sinTheta,Ene);
+
+  //  G4double theta = u * electron_mass_c2 / photonEnergy ;
+  //  G4double phi  = twopi * G4UniformRand() ;
+
+  G4double phi,psi=0.;
+
+  //corrected e+ e- angular angular distribution //preliminary!
+
+  //  if(photonEnergy>50*MeV)
+  // {
+  phi = SetPhi(photonEnergy);
+  psi = SetPsi(photonEnergy,phi);
+  //  }
+  //else
+  // {
+  //psi = G4UniformRand()*2.*pi;
+  //phi = pi; // coplanar
+  // }
+
+  Psi = psi;
+  Phi = phi;
+  //G4cout << "PHI " << phi << G4endl;
+  //G4cout << "PSI " << psi << G4endl;
+
+  G4double phie = psi; //azimuthal angle for the electron
+
+  G4double dirX = sinTheta*cos(phie);
+  G4double dirY = sinTheta*sin(phie);
+  G4double dirZ = cosTheta;
+  G4ThreeVector electronDirection(dirX,dirY,dirZ);
+  // Kinematics of the created pair:
+  // the electron and positron are assumed to have a symetric angular
+  // distribution with respect to the Z axis along the parent photon
+
+  //G4double localEnergyDeposit = 0. ;
+
+  G4double electronKineEnergy = std::max(0.,electronTotEnergy - electron_mass_c2) ;
+
+  SystemOfRefChange(gammaDirection0,electronDirection,
+		    gammaPolarization0);
+
+  G4DynamicParticle* particle1 = new G4DynamicParticle (G4Electron::Electron(),
+							electronDirection,
+							electronKineEnergy);
+
+  // The e+ is always created (even with kinetic energy = 0) for further annihilation
+
+  Ene = positronTotEnergy/electron_mass_c2; // Normalized energy
+
+  cosTheta = 0.;
+  sinTheta = 0.;
+
+  SetTheta(&cosTheta,&sinTheta,Ene);
+  G4double phip = phie+phi; //azimuthal angle for the positron
+
+  dirX = sinTheta*cos(phip);
+  dirY = sinTheta*sin(phip);
+  dirZ = cosTheta;
+  G4ThreeVector positronDirection(dirX,dirY,dirZ);
+
+  G4double positronKineEnergy = std::max(0.,positronTotEnergy - electron_mass_c2) ;
+  SystemOfRefChange(gammaDirection0,positronDirection,
+		    gammaPolarization0);
+
+  // Create G4DynamicParticle object for the particle2
+  G4DynamicParticle* particle2 = new G4DynamicParticle(G4Positron::Positron(),
+                                                       positronDirection, positronKineEnergy);
+
+
+  fvect->push_back(particle1);
+  fvect->push_back(particle2);
+
+
+
+  // Kill the incident photon
+
+
+
+  // Create lists of pointers to DynamicParticles (photons and electrons)
+  // (Is the electron vector necessary? To be checked)
+  //  std::vector<G4DynamicParticle*>* photonVector = 0;
+  //std::vector<G4DynamicParticle*> electronVector;
+
+  fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
+  fParticleChange->SetProposedKineticEnergy(0.);
+  fParticleChange->ProposeTrackStatus(fStopAndKill);
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction1(G4double screenVariable)
+{
+  // Compute the value of the screening function 3*phi1 - phi2
+
+  G4double value;
+
+  if (screenVariable > 1.)
+    value = 42.24 - 8.368 * log(screenVariable + 0.952);
+  else
+    value = 42.392 - screenVariable * (7.796 - 1.961 * screenVariable);
+
+  return value;
+}
+
+
+
+G4double G4LivermorePolarizedGammaConversionModel::ScreenFunction2(G4double screenVariable)
+{
+  // Compute the value of the screening function 1.5*phi1 - 0.5*phi2
+
+  G4double value;
+
+  if (screenVariable > 1.)
+    value = 42.24 - 8.368 * log(screenVariable + 0.952);
+  else
+    value = 41.405 - screenVariable * (5.828 - 0.8945 * screenVariable);
+
+  return value;
+}
+
+
+void G4LivermorePolarizedGammaConversionModel::SetTheta(G4double* p_cosTheta, G4double* p_sinTheta, G4double Energy)
+{
+
+  // to avoid computational errors since Theta could be very small
+  // Energy in Normalized Units (!)
+
+  G4double Momentum = sqrt(Energy*Energy -1);
+  G4double Rand = G4UniformRand();
+
+  *p_cosTheta = (Energy*((2*Rand)- 1) + Momentum)/((Momentum*(2*Rand-1))+Energy);
+  *p_sinTheta = (2*sqrt(Rand*(1-Rand)))/(Momentum*(2*Rand-1)+Energy);
+}
+
+
+
+G4double G4LivermorePolarizedGammaConversionModel::SetPhi(G4double Energy)
+{
+
+
+  G4double value = 0.;
+  G4double Ene = Energy/MeV;
+
+  G4double pl[4];
+
+
+  G4double pt[2];
+  G4double xi = 0;
+  G4double xe = 0.;
+  G4double n1=0.;
+  G4double n2=0.;
+
+
+  if (Ene>=50.)
+    {
+      const G4double ay0=5.6, by0=18.6, aa0=2.9, ba0 = 8.16E-3;
+      const G4double aw = 0.0151, bw = 10.7, cw = -410.;
+
+      const G4double axc = 3.1455, bxc = -1.11, cxc = 310.;
+
+      pl[0] = Fln(ay0,by0,Ene);
+      pl[1] = aa0 + ba0*(Ene);
+      pl[2] = Poli(aw,bw,cw,Ene);
+      pl[3] = Poli(axc,bxc,cxc,Ene);
+
+      const G4double abf = 3.1216, bbf = 2.68;
+      pt[0] = -1.4;
+      pt[1] = abf + bbf/Ene;
+
+
+
+      //G4cout << "PL > 50. "<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
+
+      xi = 3.0;
+      xe = Encu(pl,pt,xi);
+      //G4cout << "ENCU "<< xe << G4endl;
+      n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
+      n2 = Finttan(pt,xe) - Finttan(pt,0.);
+    }
+  else
+    {
+      const G4double ay0=0.144, by0=0.11;
+      const G4double aa0=2.7, ba0 = 2.74;
+      const G4double aw = 0.21, bw = 10.8, cw = -58.;
+      const G4double axc = 3.17, bxc = -0.87, cxc = -6.;
+
+      pl[0] = Fln(ay0, by0, Ene);
+      pl[1] = Fln(aa0, ba0, Ene);
+      pl[2] = Poli(aw,bw,cw,Ene);
+      pl[3] = Poli(axc,bxc,cxc,Ene);
+
+      //G4cout << "PL < 50."<< pl[0] << " " << pl[1] << " " << pl[2] << " " <<pl[3] << " " << G4endl;
+      //G4cout << "ENCU "<< xe << G4endl;
+      n1 = Fintlor(pl,pi) - Fintlor(pl,xe);
+
+    }
+
+
+  G4double n=0.;
+  n = n1+n2;
+
+  G4double c1 = 0.;
+  c1 = Glor(pl, xe);
+
+  G4double xm = 0.;
+  xm = Flor(pl,pl[3])*Glor(pl,pl[3]);
+
+  G4double r1,r2,r3;
+  G4double xco=0.;
+
+  if (Ene>=50.)
+    {
+      r1= G4UniformRand();
+      if( r1>=n2/n)
+        {
+          do
+	    {
+              r2 = G4UniformRand();
+              value = Finvlor(pl,xe,r2);
+              xco = Glor(pl,value)/c1;
+              r3 = G4UniformRand();
+            } while(r3>=xco);
+        }
+      else
+        {
+          value = Finvtan(pt,n,r1);
+        }
+    }
+  else
+    {
+      do
+        {
+          r2 = G4UniformRand();
+          value = Finvlor(pl,xe,r2);
+          xco = Glor(pl,value)/c1;
+          r3 = G4UniformRand();
+        } while(r3>=xco);
+    }
+
+  //  G4cout << "PHI = " <<value <<  G4endl;
+  return value;
+}
+G4double G4LivermorePolarizedGammaConversionModel::SetPsi(G4double Energy, G4double Phi)
+{
+
+  G4double value = 0.;
+  G4double Ene = Energy/MeV;
+
+  G4double p0l[4];
+  G4double ppml[4];
+  G4double p0t[2];
+  G4double ppmt[2];
+
+  G4double xi = 0.;
+  G4double xe0 = 0.;
+  G4double xepm = 0.;
+
+  if (Ene>=50.)
+    {
+      const G4double ay00 = 3.4, by00 = 9.8, aa00 = 1.34, ba00 = 5.3;
+      const G4double aw0 = 0.014, bw0 = 9.7, cw0 = -2.E4;
+      const G4double axc0 = 3.1423, bxc0 = -2.35, cxc0 = 0.;
+      const G4double ay0p = 1.53, by0p = 3.2, aap = 0.67, bap = 8.5E-3;
+      const G4double awp = 6.9E-3, bwp = 12.6, cwp = -3.8E4;
+      const G4double axcp = 2.8E-3,bxcp = -3.133;
+      const G4double abf0 = 3.1213, bbf0 = 2.61;
+      const G4double abfpm = 3.1231, bbfpm = 2.84;
+
+      p0l[0] = Fln(ay00, by00, Ene);
+      p0l[1] = Fln(aa00, ba00, Ene);
+      p0l[2] = Poli(aw0, bw0, cw0, Ene);
+      p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
+
+      ppml[0] = Fln(ay0p, by0p, Ene);
+      ppml[1] = aap + bap*(Ene);
+      ppml[2] = Poli(awp, bwp, cwp, Ene);
+      ppml[3] = Fln(axcp,bxcp,Ene);
+
+      p0t[0] = -0.81;
+      p0t[1] = abf0 + bbf0/Ene;
+      ppmt[0] = -0.6;
+      ppmt[1] = abfpm + bbfpm/Ene;
+
+      //G4cout << "P0L > 50"<< p0l[0] << " " << p0l[1] << " " << p0l[2] << " " <<p0l[3] << " " << G4endl;
+      //G4cout << "PPML > 50"<< ppml[0] << " " << ppml[1] << " " << ppml[2] << " " <<ppml[3] << " " << G4endl;
+
+      xi = 3.0;
+      xe0 = Encu(p0l, p0t, xi);
+      //G4cout << "ENCU1 "<< xe0 << G4endl;
+      xepm = Encu(ppml, ppmt, xi);
+
+
+      //G4cout << "ENCU2 "<< xepm << G4endl;
+
+    }
+  else
+    {
+      const G4double ay00 = 2.82, by00 = 6.35;
+      const G4double aa00 = -1.75, ba00 = 0.25;
+
+      const G4double aw0 = 0.028, bw0 = 5., cw0 = -50.;
+      const G4double axc0 = 3.14213, bxc0 = -2.3, cxc0 = 5.7;
+      const G4double ay0p = 1.56, by0p = 3.6;
+      const G4double aap = 0.86, bap = 8.3E-3;
+      const G4double awp = 0.022, bwp = 7.4, cwp = -51.;
+      const G4double xcp = 3.1486;
+
+      p0l[0] = Fln(ay00, by00, Ene);
+      p0l[1] = aa00+pow(Ene, ba00);
+      p0l[2] = Poli(aw0, bw0, cw0, Ene);
+      p0l[3] = Poli(axc0, bxc0, cxc0, Ene);
+      ppml[0] = Fln(ay0p, by0p, Ene);
+      ppml[1] = aap + bap*(Ene);
+      ppml[2] = Poli(awp, bwp, cwp, Ene);
+      ppml[3] = xcp;
+
+    }
+
+
+
+  G4double a,b=0.;
+
+  if (Ene>=50.)
+    {
+      if (Phi>xepm)
+	{
+          b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,Phi));
+        }
+      else
+        {
+          b = Ftan(ppmt,Phi);
+        }
+      if (Phi>xe0)
+        {
+          a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,Phi));
+        }
+      else
+        {
+          a = Ftan(p0t,Phi);
+        }
+    }
+  else
+    {
+      b = (ppml[0]+2*ppml[1]*ppml[2]*Flor(ppml,Phi));
+      a = (p0l[0]+2*p0l[1]*p0l[2]*Flor(p0l,Phi));
+    }
+  G4double nr =0.;
+
+  if (b>a)
+    {
+      nr = 1./b;
+    }
+  else
+    {
+      nr = 1./a;
+    }
+
+  G4double r1,r2=0.;
+  G4double r3 =-1.;
+  do
+    {
+      r1 = G4UniformRand();
+      r2 = G4UniformRand();
+      value = r2*pi;
+      r3 = nr*(a*cos(value)*cos(value) + b*sin(value)*sin(value));
+    }while(r1>r3);
+
+  return value;
+}
+
+
+G4double G4LivermorePolarizedGammaConversionModel::Poli
+(G4double a, G4double b, G4double c, G4double x)
+{
+  G4double value=0.;
+  if(x>0.)
+    {
+      value =(a + b/x + c/(x*x*x));
+    }
+  else
+    {
+      //G4cout << "ERROR in Poli! " << G4endl;
+    }
+  return value;
+}
+G4double G4LivermorePolarizedGammaConversionModel::Fln
+(G4double a, G4double b, G4double x)
+{
+  G4double value=0.;
+  if(x>0.)
+    {
+      value =(a*log(x)-b);
+    }
+  else
+    {
+      //G4cout << "ERROR in Fln! " << G4endl;
+    }
+  return value;
+}
+
+
+G4double G4LivermorePolarizedGammaConversionModel::Encu
+(G4double* p_p1, G4double* p_p2, G4double x)
+{
+  G4double value=0.;
+  G4int i=0;
+  G4double fx = 1.;
+
+  do
+    {
+      x -= (Flor(p_p1, x)*Glor(p_p1,x) - Ftan(p_p2, x))/
+        (Fdlor(p_p1,x) - Fdtan(p_p2,x));
+      fx = Flor(p_p1,x)*Glor(p_p1,x) - Ftan(p_p2, x);
+      i += 1;
+      //G4cout << abs(fx) << " " << i << " " << x << "dentro ENCU " << G4endl;
+    } while( (i<100) && (abs(fx) > 1e-6)) ;
+
+  if (i>100||x>pi) x = 3.0;
+  value = x;
+
+  if (value<0.) value=0.;
+
+  return value;
+}
+
+
+
+
+G4double G4LivermorePolarizedGammaConversionModel::Flor(G4double* p_p1, G4double x)
+{
+  G4double value =0.;
+  // G4double y0 = p_p1[0];
+  // G4double A = p_p1[1];
+  G4double w = p_p1[2];
+  G4double xc = p_p1[3];
+
+  value = 1./(pi*(w*w + 4.*(x-xc)*(x-xc)));
+  return value;
+}
+
+G4double G4LivermorePolarizedGammaConversionModel::Glor(G4double* p_p1, G4double x)
+{
+  G4double value =0.;
+  G4double y0 = p_p1[0];
+  G4double A = p_p1[1];
+  G4double w = p_p1[2];
+  G4double xc = p_p1[3];
+
+  value = (y0 *pi*(w*w +  4.*(x-xc)*(x-xc)) + 2.*A*w);
+  return value;
+}
+
+G4double G4LivermorePolarizedGammaConversionModel::Fdlor(G4double* p_p1, G4double x)
+{
+  G4double value =0.;
+  //G4double y0 = p_p1[0];
+  G4double A = p_p1[1];
+  G4double w = p_p1[2];
+  G4double xc = p_p1[3];
+
+  value = (-16.*A*w*(x-xc))/
+    (pi*(w*w+4.*(x-xc)*(x-xc))*(w*w+4.*(x-xc)*(x-xc)));
+  return value;
+}
+
+
+G4double G4LivermorePolarizedGammaConversionModel::Fintlor(G4double* p_p1, G4double x)
+{
+  G4double value =0.;
+  G4double y0 = p_p1[0];
+  G4double A = p_p1[1];
+  G4double w = p_p1[2];
+  G4double xc = p_p1[3];
+
+  value = y0*x + A*atan( 2*(x-xc)/w) / pi;
+  return value;
+}
+
+
+G4double G4LivermorePolarizedGammaConversionModel::Finvlor(G4double* p_p1, G4double x, G4double r)
+{
+  G4double value = 0.;
+  G4double nor = 0.;
+  //G4double y0 = p_p1[0];
+  //  G4double A = p_p1[1];
+  G4double w = p_p1[2];
+  G4double xc = p_p1[3];
+
+  nor = atan(2.*(pi-xc)/w)/(2.*pi*w) - atan(2.*(x-xc)/w)/(2.*pi*w);
+  value = xc - (w/2.)*tan(-2.*r*nor*pi*w+atan(2*(xc-x)/w));
+
+  return value;
+}
+
+G4double G4LivermorePolarizedGammaConversionModel::Ftan(G4double* p_p1, G4double x)
+{
+  G4double value =0.;
+  G4double a = p_p1[0];
+  G4double b = p_p1[1];
+
+  value = a /(x-b);
+  return value;
+}
+
+
+G4double G4LivermorePolarizedGammaConversionModel::Fdtan(G4double* p_p1, G4double x)
+{
+  G4double value =0.;
+  G4double a = p_p1[0];
+  G4double b = p_p1[1];
+
+  value = -1.*a / ((x-b)*(x-b));
+  return value;
+}
+
+
+G4double G4LivermorePolarizedGammaConversionModel::Finttan(G4double* p_p1, G4double x)
+{
+  G4double value =0.;
+  G4double a = p_p1[0];
+  G4double b = p_p1[1];
+
+
+  value = a*log(b-x);
+  return value;
+}
+
+G4double G4LivermorePolarizedGammaConversionModel::Finvtan(G4double* p_p1, G4double cnor, G4double r)
+{
+  G4double value =0.;
+  G4double a = p_p1[0];
+  G4double b = p_p1[1];
+
+  value = b*(1-exp(r*cnor/a));
+
+  return value;
+}
+
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4ThreeVector G4LivermorePolarizedGammaConversionModel::SetPerpendicularVector(G4ThreeVector& a)
+{
+  G4double dx = a.x();
+  G4double dy = a.y();
+  G4double dz = a.z();
+  G4double x = dx < 0.0 ? -dx : dx;
+  G4double y = dy < 0.0 ? -dy : dy;
+  G4double z = dz < 0.0 ? -dz : dz;
+  if (x < y) {
+    return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
+  }else{
+    return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetRandomPolarization(G4ThreeVector& direction0)
+{
+  G4ThreeVector d0 = direction0.unit();
+  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
+  G4ThreeVector a0 = a1.unit(); // unit vector
+
+  G4double rand1 = G4UniformRand();
+  
+  G4double angle = twopi*rand1; // random polar angle
+  G4ThreeVector b0 = d0.cross(a0); // cross product
+  
+  G4ThreeVector c;
+  
+  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
+  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
+  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
+  
+  G4ThreeVector c0 = c.unit();
+
+  return c0;
+  
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4ThreeVector G4LivermorePolarizedGammaConversionModel::GetPerpendicularPolarization
+(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
+{
+
+  // 
+  // The polarization of a photon is always perpendicular to its momentum direction.
+  // Therefore this function removes those vector component of gammaPolarization, which
+  // points in direction of gammaDirection
+  //
+  // Mathematically we search the projection of the vector a on the plane E, where n is the
+  // plains normal vector.
+  // The basic equation can be found in each geometry book (e.g. Bronstein):
+  // p = a - (a o n)/(n o n)*n
+  
+  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+
+void G4LivermorePolarizedGammaConversionModel::SystemOfRefChange
+    (G4ThreeVector& direction0,G4ThreeVector& direction1,
+     G4ThreeVector& polarization0)
+{
+  // direction0 is the original photon direction ---> z
+  // polarization0 is the original photon polarization ---> x
+  // need to specify y axis in the real reference frame ---> y 
+  G4ThreeVector Axis_Z0 = direction0.unit();
+  G4ThreeVector Axis_X0 = polarization0.unit();
+  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
+  
+  G4double direction_x = direction1.getX();
+  G4double direction_y = direction1.getY();
+  G4double direction_z = direction1.getZ();
+  
+  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 +  direction_z*Axis_Z0).unit();
+  
+}
+
+
+
+
Index: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedPhotoElectricModel.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedPhotoElectricModel.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedPhotoElectricModel.cc	(revision 1197)
@@ -0,0 +1,553 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+
+#include "G4LivermorePolarizedPhotoElectricModel.hh"
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+using namespace std;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4LivermorePolarizedPhotoElectricModel::G4LivermorePolarizedPhotoElectricModel(const G4ParticleDefinition*,
+                                             const G4String& nam)
+  :G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),crossSectionHandler(0), shellCrossSectionHandler(0)
+{
+  lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ?
+  highEnergyLimit = 100 * GeV;
+  SetLowEnergyLimit(lowEnergyLimit);
+  SetHighEnergyLimit(highEnergyLimit);
+  
+  verboseLevel= 0;
+  // Verbosity scale:
+  // 0 = nothing 
+  // 1 = warning for energy non-conservation 
+  // 2 = details of energy budget
+  // 3 = calculation of cross sections, file openings, sampling of atoms
+  // 4 = entering in methods
+
+  G4cout << "Livermore Polarized PhotoElectric is constructed " << G4endl
+         << "Energy range: "
+         << lowEnergyLimit / keV << " keV - "
+         << highEnergyLimit / GeV << " GeV"
+         << G4endl;
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4LivermorePolarizedPhotoElectricModel::~G4LivermorePolarizedPhotoElectricModel()
+{  
+  if (meanFreePathTable)   delete meanFreePathTable;
+  if (crossSectionHandler) delete crossSectionHandler;
+  if (shellCrossSectionHandler) delete shellCrossSectionHandler;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LivermorePolarizedPhotoElectricModel::Initialise(const G4ParticleDefinition* particle,
+                                       const G4DataVector& cuts)
+{
+  if (verboseLevel > 3)
+    G4cout << "Calling G4LivermorePolarizedPhotoElectricModel::Initialise()" << G4endl;
+
+  if (crossSectionHandler)
+  {
+    crossSectionHandler->Clear();
+    delete crossSectionHandler;
+  }
+
+  if (shellCrossSectionHandler)
+    {
+      shellCrossSectionHandler->Clear();
+      delete shellCrossSectionHandler;
+    }
+
+
+  // Energy limits
+  
+  if (LowEnergyLimit() < lowEnergyLimit)
+    {
+      G4cout << "G4LivermorePolarizedPhotoElectricModel: low energy limit increased from " << 
+	LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
+      SetLowEnergyLimit(lowEnergyLimit);
+    }
+
+  if (HighEnergyLimit() > highEnergyLimit)
+    {
+      G4cout << "G4LivermorePolarizedPhotoElectricModel: high energy limit decreased from " << 
+	HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
+      SetHighEnergyLimit(highEnergyLimit);
+    }
+
+  // Reading of data files - all materials are read
+  
+  crossSectionHandler = new G4CrossSectionHandler;
+  crossSectionHandler->Clear();
+  G4String crossSectionFile = "phot/pe-cs-";
+  crossSectionHandler->LoadData(crossSectionFile);
+
+  meanFreePathTable = 0;
+  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
+
+  shellCrossSectionHandler = new G4CrossSectionHandler();
+  shellCrossSectionHandler->Clear();
+  G4String shellCrossSectionFile = "phot/pe-ss-cs-";
+  shellCrossSectionHandler->LoadShellData(shellCrossSectionFile);
+
+
+  //
+  if (verboseLevel > 2) 
+    G4cout << "Loaded cross section files for Livermore Polarized PhotoElectric model" << G4endl;
+
+  InitialiseElementSelectors(particle,cuts);
+
+  G4cout << "Livermore Polarized PhotoElectric model is initialized " << G4endl
+         << "Energy range: "
+         << LowEnergyLimit() / keV << " keV - "
+         << HighEnergyLimit() / GeV << " GeV"
+         << G4endl;
+
+  //
+    
+  if(isInitialised) return;
+
+  if(pParticleChange)
+    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
+  else
+    fParticleChange = new G4ParticleChangeForGamma();
+    
+  isInitialised = true;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4LivermorePolarizedPhotoElectricModel::ComputeCrossSectionPerAtom(
+                                       const G4ParticleDefinition*,
+                                             G4double GammaEnergy,
+                                             G4double Z, G4double,
+                                             G4double, G4double)
+{
+  if (verboseLevel > 3)
+    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedPhotoElectricModel" << G4endl;
+
+  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
+  return cs;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LivermorePolarizedPhotoElectricModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
+							       const G4MaterialCutsCouple* couple,
+							       const G4DynamicParticle* aDynamicGamma,
+							       G4double,
+							       G4double)
+{
+
+  // Fluorescence generated according to:
+  // J. Stepanek ,"A program to determine the radiation spectra due to a single atomic
+  // subshell ionisation by a particle or due to deexcitation or decay of radionuclides",
+  // Comp. Phys. Comm. 1206 pp 1-1-9 (1997)
+
+  if (verboseLevel > 3)
+    G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedPhotoElectricModel" << G4endl;
+
+  G4double photonEnergy = aDynamicGamma->GetKineticEnergy();
+  // Within energy limit?
+
+  if(photonEnergy <= lowEnergyLimit)
+    {
+      fParticleChange->ProposeTrackStatus(fStopAndKill);
+      fParticleChange->SetProposedKineticEnergy(0.);
+      fParticleChange->ProposeLocalEnergyDeposit(photonEnergy);
+      return;
+    }
+
+
+  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
+
+  // Protection: a polarisation parallel to the
+  // direction causes problems;
+  // in that case find a random polarization
+
+  G4ThreeVector photonDirection = aDynamicGamma->GetMomentumDirection();
+
+  // Make sure that the polarization vector is perpendicular to the
+  // gamma direction. If not
+
+  if(!(gammaPolarization0.isOrthogonal(photonDirection, 1e-6))||(gammaPolarization0.mag()==0))
+    { // only for testing now
+      gammaPolarization0 = GetRandomPolarization(photonDirection);
+    }
+  else
+    {
+      if ( gammaPolarization0.howOrthogonal(photonDirection) != 0)
+	{
+	  gammaPolarization0 = GetPerpendicularPolarization(photonDirection, gammaPolarization0);
+	}
+    }
+
+  // End of Protection
+
+  //  G4double E0_m = photonEnergy / electron_mass_c2 ;
+
+  // Select randomly one element in the current material
+
+  G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy);
+
+  // Select the ionised shell in the current atom according to shell cross sections
+
+  size_t shellIndex = shellCrossSectionHandler->SelectRandomShell(Z,photonEnergy);
+
+  const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
+  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
+  G4double bindingEnergy = shell->BindingEnergy();
+  G4int shellId = shell->ShellId();
+
+  // Create lists of pointers to DynamicParticles (photons and electrons)
+  // (Is the electron vector necessary? To be checked)
+  std::vector<G4DynamicParticle*>* photonVector = 0;
+  std::vector<G4DynamicParticle*> electronVector;
+
+  G4double energyDeposit = 0.0;
+
+  // Primary outgoing  electron
+
+  G4double eKineticEnergy = photonEnergy - bindingEnergy;
+
+
+  if (eKineticEnergy > 0.)
+    { 
+
+      G4double costheta = SetCosTheta(eKineticEnergy);
+      G4double sintheta = sqrt(1. - costheta*costheta);
+      G4double phi = SetPhi(photonEnergy,eKineticEnergy,costheta);
+      G4double dirX = sintheta*cos(phi);
+      G4double dirY = sintheta*sin(phi);
+      G4double dirZ = costheta;
+      G4ThreeVector electronDirection(dirX, dirY, dirZ);
+      SystemOfRefChange(photonDirection, electronDirection, gammaPolarization0);
+      G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(), 
+							   electronDirection, 
+							   eKineticEnergy);
+      electronVector.push_back(electron);
+    }
+  else
+    {
+      bindingEnergy = photonEnergy;
+    }
+  
+
+  G4int nElectrons = electronVector.size();
+  size_t nTotPhotons = 0;
+  G4int nPhotons=0;
+
+  const G4ProductionCutsTable* theCoupleTable=
+    G4ProductionCutsTable::GetProductionCutsTable();
+  size_t index = couple->GetIndex();
+  G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[index];
+  cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
+
+  G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[index];
+  cute = std::min(cutForLowEnergySecondaryPhotons,cute);
+
+  G4DynamicParticle* aPhoton;
+
+  // Generation of fluorescence
+  // Data in EADL are available only for Z > 5
+  // Protection to avoid generating photons in the unphysical case of
+  // shell binding energy > photon energy
+  if (Z > 5  && (bindingEnergy > cutg || bindingEnergy > cute))
+    {
+      photonVector = deexcitationManager.GenerateParticles(Z,shellId);
+      nTotPhotons = photonVector->size();
+      for (size_t k=0; k<nTotPhotons; k++)
+        {
+          aPhoton = (*photonVector)[k];
+          if (aPhoton)
+            {
+              G4double itsCut = cutg;
+              if(aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
+
+	      G4double itsEnergy = aPhoton->GetKineticEnergy();
+
+              if (itsEnergy > itsCut && itsEnergy <= bindingEnergy)
+                {
+                  nPhotons++;
+                  // Local energy deposit is given as the sum of the
+                  // energies of incident photons minus the energies
+                  // of the outcoming fluorescence photons
+                  bindingEnergy -= itsEnergy;
+
+                }
+	      else
+                {
+                  delete aPhoton;
+                  (*photonVector)[k] = 0;
+                }
+	    }
+	} 
+    }      
+  energyDeposit += bindingEnergy;
+  // Final state
+
+  for (G4int l = 0; l<nElectrons; l++ )
+    {
+      aPhoton = electronVector[l];
+      if(aPhoton) {
+        fvect->push_back(aPhoton);
+      }
+    }
+  for ( size_t ll = 0; ll < nTotPhotons; ll++)
+    {
+      aPhoton = (*photonVector)[ll];
+      if(aPhoton) {
+        fvect->push_back(aPhoton);
+      }
+    }
+
+  delete photonVector;
+
+  if (energyDeposit < 0)
+    {
+      G4cout << "WARNING - "
+             << "G4LowEnergyPhotoElectric::PostStepDoIt - Negative energy deposit"
+             << G4endl;
+      energyDeposit = 0;
+    }
+
+  // kill incident photon
+  fParticleChange->ProposeMomentumDirection( 0., 0., 0. );
+  fParticleChange->SetProposedKineticEnergy(0.);
+  fParticleChange->ProposeTrackStatus(fStopAndKill);
+  fParticleChange->ProposeLocalEnergyDeposit(energyDeposit);
+
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LivermorePolarizedPhotoElectricModel::SetCutForLowEnSecPhotons(G4double cut)
+{
+  cutForLowEnergySecondaryPhotons = cut;
+  deexcitationManager.SetCutForSecondaryPhotons(cut);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LivermorePolarizedPhotoElectricModel::SetCutForLowEnSecElectrons(G4double cut)
+{
+  cutForLowEnergySecondaryElectrons = cut;
+  deexcitationManager.SetCutForAugerElectrons(cut);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4LivermorePolarizedPhotoElectricModel::ActivateAuger(G4bool val)
+{
+  deexcitationManager.ActivateAugerElectronProduction(val);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+
+G4double G4LivermorePolarizedPhotoElectricModel::SetCosTheta(G4double energyE)
+
+{
+  G4double rand1,rand2,onemcost,greject;
+  G4double masarep = 510.99906*keV;
+
+  G4double gamma = 1. + energyE/masarep;
+  G4double gamma2 = gamma*gamma;
+
+  G4double beta = sqrt((gamma2 - 1.)/gamma2);
+
+  G4double alfa = 1./beta - 1.;
+
+  G4double g1 = 0.5*beta*gamma*(gamma-1.)*(gamma-2.);
+
+  G4double alfap2 = alfa+2.;
+
+  G4double grejectmax = 2.*(g1+1./alfa);
+
+  do
+    {
+      rand1 = G4UniformRand();
+      onemcost = 2.*alfa*(2.*rand1 + alfap2 * sqrt(rand1))/
+	(alfap2*alfap2 - 4.*rand1);
+      greject = (2. - onemcost)*(g1+1./(alfa+onemcost));
+      rand2 = G4UniformRand();
+    }
+  while (rand2*grejectmax > greject);
+  G4double cosTheta = 1. - onemcost;
+  return cosTheta;
+} 
+
+
+
+G4double G4LivermorePolarizedPhotoElectricModel::SetPhi(G4double Ph_energy,
+						   G4double E_energy,
+						   G4double costheta)
+{
+  G4double epsilon = E_energy/electron_mass_c2;
+  G4double k = Ph_energy/electron_mass_c2;
+  G4double gamma = 1. + epsilon;
+  G4double gamma2 = gamma*gamma;
+  G4double beta = sqrt((gamma2 - 1.)/gamma2);
+
+  G4double d = (2./(k*gamma*(1-beta*costheta))-1)*(1/k);
+
+  G4double norm_factor = 1 +2*d;
+
+  G4double rnd1; 
+  G4double rnd2;
+  G4double phi, phiprob;
+
+  do
+    {
+      rnd1 =G4UniformRand();
+      rnd2 =G4UniformRand();
+      phi = rnd1*twopi;
+      phiprob = 1 +2*d*cos(phi)*cos(phi);
+    }
+  while (rnd2*norm_factor > phiprob);
+  return phi;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4ThreeVector G4LivermorePolarizedPhotoElectricModel::SetPerpendicularVector(G4ThreeVector& a)
+{
+  G4double dx = a.x();
+  G4double dy = a.y();
+  G4double dz = a.z();
+  G4double x = dx < 0.0 ? -dx : dx;
+  G4double y = dy < 0.0 ? -dy : dy;
+  G4double z = dz < 0.0 ? -dz : dz;
+  if (x < y) {
+    return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
+  }else{
+    return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4ThreeVector G4LivermorePolarizedPhotoElectricModel::GetRandomPolarization(G4ThreeVector& direction0)
+{
+  G4ThreeVector d0 = direction0.unit();
+  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
+  G4ThreeVector a0 = a1.unit(); // unit vector
+
+  G4double rand1 = G4UniformRand();
+  
+  G4double angle = twopi*rand1; // random polar angle
+  G4ThreeVector b0 = d0.cross(a0); // cross product
+  
+  G4ThreeVector c;
+  
+  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
+  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
+  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
+  
+  G4ThreeVector c0 = c.unit();
+
+  return c0;
+  
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4ThreeVector G4LivermorePolarizedPhotoElectricModel::GetPerpendicularPolarization
+(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
+{
+
+  // 
+  // The polarization of a photon is always perpendicular to its momentum direction.
+  // Therefore this function removes those vector component of gammaPolarization, which
+  // points in direction of gammaDirection
+  //
+  // Mathematically we search the projection of the vector a on the plane E, where n is the
+  // plains normal vector.
+  // The basic equation can be found in each geometry book (e.g. Bronstein):
+  // p = a - (a o n)/(n o n)*n
+  
+  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+
+void G4LivermorePolarizedPhotoElectricModel::SystemOfRefChange
+    (G4ThreeVector& direction0,G4ThreeVector& direction1,
+     G4ThreeVector& polarization0)
+{
+  // direction0 is the original photon direction ---> z
+  // polarization0 is the original photon polarization ---> x
+  // need to specify y axis in the real reference frame ---> y 
+  G4ThreeVector Axis_Z0 = direction0.unit();
+  G4ThreeVector Axis_X0 = polarization0.unit();
+  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
+  
+  G4double direction_x = direction1.getX();
+  G4double direction_y = direction1.getY();
+  G4double direction_z = direction1.getZ();
+  
+  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 +  direction_z*Axis_Z0).unit();
+  
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4LivermorePolarizedPhotoElectricModel::GetMeanFreePath(const G4Track& track,
+						      G4double,
+						      G4ForceCondition*)
+{
+
+  const G4DynamicParticle* photon = track.GetDynamicParticle();
+  G4double energy = photon->GetKineticEnergy();
+  G4Material* material = track.GetMaterial();
+  //  size_t materialIndex = material->GetIndex();
+
+  G4double meanFreePath = DBL_MAX;
+
+  //  if (energy > highEnergyLimit)
+  //    meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
+  //  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
+  //  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
+
+  G4double cross = shellCrossSectionHandler->ValueForMaterial(material,energy);
+  if(cross > 0.0) meanFreePath = 1.0/cross;
+
+  return meanFreePath;
+
+
+}
+
+
Index: trunk/source/processes/electromagnetic/lowenergy/src/G4OrlicLiCrossSection.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4OrlicLiCrossSection.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4OrlicLiCrossSection.cc	(revision 1197)
@@ -0,0 +1,469 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//$Id: G4OrlicLiCrossSection.cc,v 1.4 2009/11/11 09:14:53 mantero Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// Author: Haifa Ben Abdelouahed
+//         
+//
+// History:
+// -----------
+//  23 Apr 2008   H. Ben Abdelouahed   1st implementation
+//  28 Apr 2008   MGP        Major revision according to a design iteration
+//  21 Apr 2009	  ALF Some correction for compatibility to G4VShellCrossSection
+//		  and changed name to G4OrlicLiCrossSection 
+//  11 Nov 2009	  ALF update and code cleaning for the Dec Release
+//
+// -------------------------------------------------------------------
+// Class description:
+// Low Energy Electromagnetic Physics, Cross section, proton ionisation, L shell
+// Further documentation available from http://www.ge.infn.it/geant4/lowE
+// -------------------------------------------------------------------
+
+
+#include "globals.hh"
+#include "G4OrlicLiCrossSection.hh"
+#include "G4Proton.hh"
+
+
+G4OrlicLiCrossSection::G4OrlicLiCrossSection()
+{ 
+
+  transitionManager =  G4AtomicTransitionManager::Instance();
+
+}
+
+G4OrlicLiCrossSection::~G4OrlicLiCrossSection()
+{ 
+
+}
+
+//this L-CrossSection calculation method is done according to 
+//I.ORLIC, C.H.SOW and S.M.TANG,International Journal of PIXE.Vol.4(1997) 217-230	
+ 
+
+//*****************************************************************************************************************************************
+
+G4double G4OrlicLiCrossSection::CalculateL1CrossSection(G4int zTarget, G4double energyIncident)
+  			       
+{
+
+  if ( (energyIncident < 0.1*MeV) || energyIncident > 10*MeV )
+
+    {return 0;}
+
+
+
+  G4double  massIncident; 
+
+  G4Proton* aProtone = G4Proton::Proton();
+    
+   massIncident = aProtone->GetPDGMass(); 
+ 
+  G4double l1BindingEnergy = (transitionManager->Shell(zTarget,1)->BindingEnergy())/keV;
+ 
+  G4double lamda =  massIncident/electron_mass_c2;
+
+  G4double normalizedEnergy =  (energyIncident/keV)/(lamda*l1BindingEnergy);
+
+  G4double x = std::log(normalizedEnergy);
+
+  G4double a0 = 0.;
+  G4double a1 = 0.;
+  G4double a2 = 0.;
+  G4double a3 = 0.;
+  G4double a4 = 0.;
+  G4double a5 = 0.;
+  G4double a6 = 0.;  
+  G4double a7 = 0.;  
+  G4double a8 = 0.; 
+  G4double a9 = 0.;
+
+                                       
+   if ( zTarget>=14 && zTarget<=40)  
+    {
+
+      return 0;
+      /*
+      // parameters used for calculating total L cross section
+      a0=12.5081;
+      a1=0.2177;
+      a2=-0.3758;
+      a3=0.0096;
+      a4=0.0073;
+      a5=0.0022;
+      a6=0.;
+      a7=0.;
+      a8=0.;
+      a9=0.; */
+    }
+  else 
+    { 
+       if ( zTarget>=41 &&  zTarget<=50 )
+       { 	  
+      a0=11.274881;
+      a1=-0.187401;
+      a2=-0.943341;
+      a3=-1.47817;
+      a4=-1.282343;
+      a5=-0.386544; 
+      a6=-0.037932;
+      a7=0.;
+      a8=0.;
+      a9=0.;
+	}
+     
+      else  
+	{
+	  if ( zTarget>=51 &&  zTarget<=60 ) 
+	    { 
+	      a0=11.242637;
+	      a1=-0.162515;
+	      a2=1.035774;
+	      a3=3.970908;
+	      a4=3.968233;
+	      a5=1.655714;
+	      a6=0.058885;
+	      a7=-0.155743;
+	      a8=-0.042228;
+	      a9=-0.003371; 
+	    } 
+	
+	  else  
+	    {
+	      if ( zTarget>=61 &&  zTarget<=70 ) 
+		{ 
+		  a0=6.476722;
+		  a1=-25.804787;
+		  a2=-54.061629;
+		  a3=-56.684589;
+		  a4=-33.223367;
+		  a5=-11.034979; 
+		  a6=-2.042851;
+		  a7=-0.194075;
+		  a8=-0.007252;
+		  a9=0.;
+		} 
+	      else  
+		{
+		  if ( zTarget>=71 &&  zTarget<=80 ) 
+		    { 
+		      a0=12.776794;
+		      a1=6.562907;
+		      a2=10.158703;
+		      a3=7.432592;
+		      a4=2.332036;
+		      a5=0.317946; 
+		      a6=0.014479;
+		      a7=0.;
+		      a8=0.;
+		      a9=0.;
+		    } 
+		  else  
+		    {
+		      if ( zTarget>=81 &&  zTarget<=92 ) 
+			{ 
+			  a0=28.243087;
+			  a1=50.199585;
+			  a2=58.281684;
+			  a3=34.130538;
+			  a4=10.268531;
+			  a5=1.525302; 
+			  a6=0.08835;
+			  a7=0.;
+			  a8=0.;
+			  a9=0.;
+			}
+		      else
+			{ 
+			  G4cout << "ERROR: L1 Cross-Section exist only for ZTarget between 41 and 92!!! " << G4endl;
+			
+			}
+		    }
+		}
+	    }
+	}
+      }
+      
+
+G4double analyticalFunction = a0 + (a1*x)+(a2*x*x)+(a3*std::pow(x,3))+(a4*std::pow(x,4))+(a5*std::pow(x,5))+(a6*std::pow(x,6))+
+	(a7*std::pow(x,7))+(a8*std::pow(x,8))+(a9*std::pow(x,9)); 
+
+
+
+  G4double L1crossSection =  std::exp(analyticalFunction)/(l1BindingEnergy*l1BindingEnergy);  
+
+
+  if (L1crossSection >= 0) {
+    return L1crossSection;
+  }
+  else {return 0;}
+
+}
+
+//*****************************************************************************************************************************************
+
+
+G4double G4OrlicLiCrossSection::CalculateL2CrossSection(G4int zTarget, G4double energyIncident)
+  			       
+{
+
+
+  if ( (energyIncident < 0.1*MeV) || energyIncident > 10*MeV )
+
+    {return 0;}
+
+
+  G4double  massIncident; 
+
+  G4Proton* aProtone = G4Proton::Proton();
+    
+   massIncident = aProtone->GetPDGMass(); 
+ 
+ G4double L2crossSection;
+
+ if (zTarget<41)
+   { 
+     L2crossSection =0.;
+   }
+ else
+   {
+
+ G4double l2BindingEnergy = (transitionManager->Shell(zTarget,2)->BindingEnergy())/keV;
+ 
+  G4double lamda =  massIncident/electron_mass_c2;
+
+  G4double normalizedEnergy =  (energyIncident/keV)/(lamda*l2BindingEnergy);
+
+  G4double x = std::log(normalizedEnergy);
+
+  G4double a0 = 0.;
+  G4double a1 = 0.;
+  G4double a2 = 0.;
+  G4double a3 = 0.;
+  G4double a4 = 0.;
+  G4double a5 = 0.;
+ 
+      if ( zTarget>=41 &&  zTarget<=50 ) 
+	{ 
+      a0=11.194798;
+      a1=0.178807;
+      a2=-0.449865;
+      a3=-0.063528;
+      a4=-0.015364;
+      a5=0.; 
+	}
+     
+      else  
+	{
+	  if ( zTarget>=51 &&  zTarget<=60 ) 
+	    { 
+	      a0=11.241409;
+	      a1=0.149635;
+	      a2=-0.633269;
+	      a3=-0.17834;
+	      a4=-0.034743;
+	      a5=0.006474;
+
+	    } 
+	
+	  else  
+	    {
+	      if ( zTarget>=61 &&  zTarget<=70 ) 
+		{ 
+		  a0=11.247424;
+		  a1=0.203051;
+		  a2=-0.219083;
+		  a3=0.164514;
+		  a4=0.058692;
+		  a5=0.007866; 
+		} 
+	      else  
+		{
+		  if ( zTarget>=71 &&  zTarget<=80 ) 
+		    { 
+		      a0=11.229924;
+		      a1=-0.087241;
+		      a2=-0.753908;
+		      a3=-0.181546;
+		      a4=-0.030406;
+		      a5=0.; 
+		    } 
+		  else  
+		    {
+		      if ( zTarget>=81 &&  zTarget<=92 ) 
+			{ 
+			  a0=11.586671;
+			  a1=0.730838;
+			  a2=-0.056713;
+			  a3=0.053262;
+			  a4=-0.003672;
+			  a5=0.; 
+			}
+		      else
+			{ 
+			  G4cout << "ERROR: L2 Cross-Section exist only for ZTarget between 41 and 92!!! " << G4endl;
+		       
+			}
+		    }
+		}
+	    }
+	}
+      
+
+ G4double analyticalFunction = a0 + (a1*x)+(a2*x*x)+(a3*std::pow(x,3))+(a4*std::pow(x,4))+(a5*std::pow(x,5)); 
+
+
+   L2crossSection =  std::exp(analyticalFunction)/(l2BindingEnergy*l2BindingEnergy); 
+                                                                            
+   } 
+
+  if (L2crossSection >= 0) {
+    return L2crossSection;
+  }
+  else {return 0;}
+   
+}
+
+//*****************************************************************************************************************************************
+
+
+G4double G4OrlicLiCrossSection::CalculateL3CrossSection(G4int zTarget, G4double energyIncident)
+  			       
+{
+
+  if ( (energyIncident < 0.1*MeV) || energyIncident > 10*MeV )
+
+    {return 0;}
+
+
+
+  G4double  massIncident; 
+
+  G4Proton* aProtone = G4Proton::Proton();
+    
+  massIncident = aProtone->GetPDGMass(); 
+
+ 
+ G4double L3crossSection;
+
+ if (zTarget<41)
+   { 
+     L3crossSection =0.;
+   }
+ else
+   {
+
+ G4double l3BindingEnergy = (transitionManager->Shell(zTarget,3)->BindingEnergy())/keV;
+ 
+ 
+  G4double lamda =  massIncident/electron_mass_c2;
+
+  G4double normalizedEnergy =  (energyIncident/keV)/(lamda*l3BindingEnergy);
+
+  G4double x = std::log(normalizedEnergy);
+
+
+  G4double a0 = 0.;
+  G4double a1 = 0.;
+  G4double a2 = 0.;
+  G4double a3 = 0.;
+  G4double a4 = 0.;
+ 
+      if ( zTarget>=41 &&  zTarget<=50 ) 
+	{ 
+      a0=11.91837;
+      a1=0.03064;
+      a2=-0.657644;
+      a3=-0.14532;
+      a4=-0.026059;      
+	}
+     
+      else  
+	{
+	  if ( zTarget>=51 &&  zTarget<=60 ) 
+	    { 
+	      a0=11.909485;
+	      a1=0.15918;
+	      a2=-0.588004;
+	      a3=-0.159466;
+	      a4=-0.033184;	      
+	    } 
+	
+	  else  
+	    {
+	      if ( zTarget>=61 &&  zTarget<=70 ) 
+		{ 
+		  a0=11.878472;
+		  a1=-0.137007;
+		  a2=-0.959475;
+		  a3=-0.316505;
+		  a4=-0.054154;	 
+		} 
+	      else  
+		{
+		  if ( zTarget>=71 &&  zTarget<=80 ) 
+		    { 
+		      a0=11.802538;
+		      a1=-0.371796;
+		      a2=-1.052238;
+		      a3=-0.28766;
+		      a4=-0.042608;	     
+		    } 
+		  else  
+		    {
+		      if ( zTarget>=81 &&  zTarget<=92 ) 
+			{ 
+			  a0=11.423712;
+			  a1=-1.428823;
+			  a2=-1.946979;
+			  a3=-0.585198;
+			  a4=-0.076467;			 
+			}
+		      else
+			{ 
+			  G4cout << "ERROR: L3 Cross-Section exist only for ZTarget between 41 and 92!!! " << G4endl;
+		       
+			}
+		    }
+		}
+	    }
+	}
+      
+
+ G4double analyticalFunction = a0 + (a1*x)+(a2*x*x)+(a3*std::pow(x,3))+(a4*std::pow(x,4)); 
+
+
+   L3crossSection =  std::exp(analyticalFunction)/(l3BindingEnergy*l3BindingEnergy); 
+                                                                            
+   } 
+  if (L3crossSection >= 0) {
+    return L3crossSection;
+  }
+  else {return 0;}
+
+
+}
Index: trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrKCrossSection.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrKCrossSection.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrKCrossSection.cc	(revision 1197)
@@ -0,0 +1,638 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//$Id: G4ecpssrKCrossSection.cc,v 1.7 2009/11/11 09:14:53 mantero Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// Author: Haifa Ben Abdelouahed
+//         
+//
+// History:
+// -----------
+//  21 Apr 2008   H. Ben Abdelouahed   1st implementation
+//  21 Apr 2008   MGP        Major revision according to a design iteration
+//  21 Apr 2009	  ALF Some correction for compatibility to G4VShellCrossSection
+//		  and changed name to G4ecpssrKCrossSection 
+//  11 Nov 2009   ALF update and code cleaning for the Dec Release 
+//
+// -------------------------------------------------------------------
+// Class description:
+// Low Energy Electromagnetic Physics, Cross section, p ionisation, K shell
+// Further documentation available from http://www.ge.infn.it/geant4/lowE
+// -------------------------------------------------------------------
+
+
+#include "globals.hh"
+#include "G4ecpssrKCrossSection.hh"
+#include "G4AtomicTransitionManager.hh"
+#include "G4NistManager.hh"
+#include "G4Proton.hh"
+#include "G4Alpha.hh"
+#include <math.h>
+#include <iostream>
+
+#include "G4SemiLogInterpolation.hh"
+
+
+G4ecpssrKCrossSection::G4ecpssrKCrossSection()
+{ 
+
+    // Storing FK data needed for medium velocities region
+
+    char *path = getenv("G4LEDATA");
+ 
+    if (!path)
+    G4Exception("G4ecpssrKCrossSection::CalculateCrossSection: G4LEDATA environment variable not set");
+
+    std::ostringstream fileName;
+    fileName << path << "/pixe/uf/FK.dat";
+    std::ifstream FK(fileName.str().c_str());
+
+    if (!FK) G4Exception("G4ecpssrKCrossSection::CalculateCrossSection: error opening FK data file");
+      
+    dummyVec.push_back(0.);
+
+    while(!FK.eof())
+    {
+	double x;
+	double y;
+	
+	FK>>x>>y;
+
+	//  Mandatory vector initialization
+        if (x != dummyVec.back()) 
+        { 
+          dummyVec.push_back(x); 
+          aVecMap[x].push_back(-1.);
+        }
+	  
+        FK>>FKData[x][y];
+
+        if (y != aVecMap[x].back()) aVecMap[x].push_back(y);
+          
+    }
+
+  // Storing C coefficients for high velocity formula
+
+  G4String fileC1("pixe/uf/c1");
+  tableC1 = new G4DNACrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
+  tableC1->LoadData(fileC1);
+
+  G4String fileC2("pixe/uf/c2");
+  tableC2 = new G4DNACrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
+  tableC2->LoadData(fileC2);
+
+  G4String fileC3("pixe/uf/c3");
+  tableC3 = new G4DNACrossSectionDataSet(new G4SemiLogInterpolation, 1.,1.);
+  tableC3->LoadData(fileC3);
+
+  //
+
+  verboseLevel=0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void print (G4double elem)
+{
+  G4cout << elem << " ";
+}
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4ecpssrKCrossSection::~G4ecpssrKCrossSection()
+{ 
+
+  delete tableC1;
+  delete tableC2;
+  delete tableC3;
+
+}
+
+//---------------------------------this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x)------
+
+G4double G4ecpssrKCrossSection::ExpIntFunction(G4int n,G4double x)
+
+{
+  G4int i;
+  G4int ii;
+  G4int nm1;
+  G4double a;
+  G4double b;
+  G4double c;
+  G4double d;
+  G4double del;
+  G4double fact;
+  G4double h;
+  G4double psi;
+  G4double ans = 0;
+  const G4double euler= 0.5772156649;
+  const G4int maxit= 100;
+  const G4double fpmin = 1.0e-30;
+  const G4double eps = 1.0e-7;
+  nm1=n-1;
+  if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
+  G4cout << "bad arguments in ExpIntFunction" << G4endl;
+  else {
+       if (n==0) ans=std::exp(-x)/x;
+        else {
+           if (x==0.0) ans=1.0/nm1;
+              else {
+                   if (x > 1.0) {
+                         	b=x+n;
+                        	c=1.0/fpmin;
+                        	d=1.0/b;
+				h=d;
+				for (i=1;i<=maxit;i++) {
+				  a=-i*(nm1+i);
+				  b +=2.0;
+				  d=1.0/(a*d+b);
+				  c=b+a/c;
+				  del=c*d;
+				  h *=del;
+				      if (std::fabs(del-1.0) < eps) {
+					ans=h*std::exp(-x);
+					return ans;
+				      }
+				}
+		   } else {
+		     ans = (nm1!=0 ? 1.0/nm1 : -std::log(x)-euler);
+		     fact=1.0;
+		     for (i=1;i<=maxit;i++) {
+		       fact *=-x/i;
+		       if (i !=nm1) del = -fact/(i-nm1);
+		       else {
+			 psi = -euler;
+			 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
+			 del=fact*(-std::log(x)+psi);
+		       }
+		       ans += del;
+		       if (std::fabs(del) < std::fabs(ans)*eps) return ans;
+		     }
+		   }
+	      }
+	}
+  }
+return ans;
+}
+//-----------------------------------------------------------------------------------------------------------
+
+ 
+
+G4double G4ecpssrKCrossSection::CalculateCrossSection(G4int zTarget,G4double massIncident, G4double energyIncident) 
+ 
+ //this K-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//					                
+
+{
+
+  //if (energyIncident < 150 *keV) {return 0;}
+
+  G4NistManager* massManager = G4NistManager::Instance();   
+
+  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
+
+  G4double  zIncident = 0; 
+  G4Proton* aProtone = G4Proton::Proton();
+  G4Alpha* aAlpha = G4Alpha::Alpha();
+
+  if (massIncident == aProtone->GetPDGMass() )
+  {
+
+   zIncident = (aProtone->GetPDGCharge())/eplus; 
+
+   //   G4cout << "zincident:" << zIncident << G4endl;
+    }
+  else
+    {
+      if (massIncident == aAlpha->GetPDGMass())
+	{
+	  
+	  zIncident  = (aAlpha->GetPDGCharge())/eplus; 
+	  
+	  //	  G4cout << "zincident:" << zIncident << G4endl;
+	}
+      else
+	{ 
+	  G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
+	  return 0;
+	}
+    }
+  
+  ///////////////////////////////////////////
+  //from here I will substitute this version with Seb's one. 
+  ///////////////////////////////////////////
+
+
+  if (verboseLevel>0) G4cout << "  massIncident=" << massIncident<< G4endl;
+
+  G4double kBindingEnergy = transitionManager->Shell(zTarget,0)->BindingEnergy();
+  if (verboseLevel>0) G4cout << "  kBindingEnergy=" << kBindingEnergy/eV<< G4endl;
+
+  G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
+  if (verboseLevel>0) G4cout << "  massTarget=" <<  massTarget<< G4endl;
+ 
+  G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2; //the mass of the system (projectile, target) 
+  if (verboseLevel>0) G4cout << "  systemMass=" <<  systemMass<< G4endl;
+
+  const G4double zkshell= 0.3;
+
+  G4double screenedzTarget = zTarget-zkshell; // screenedzTarget is the screened nuclear charge of the target
+
+  const G4double rydbergMeV= 13.6e-6;
+ 
+  G4double tetaK = kBindingEnergy/((screenedzTarget*screenedzTarget)*rydbergMeV);  //tetaK denotes the reduced binding energy of the electron
+  if (verboseLevel>0) G4cout << "  tetaK=" <<  tetaK<< G4endl;
+
+  G4double velocity =(2./(tetaK*screenedzTarget))*std::pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
+  if (verboseLevel>0) G4cout << "  velocity=" << velocity<< G4endl;
+
+  const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;        
+  if (verboseLevel>0) G4cout << "  bohrPow2Barn=" <<  bohrPow2Barn<< G4endl;
+
+  G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*std::pow(screenedzTarget,-4.);     //sigma0 is the initial cross section of K shell at stable state
+  if (verboseLevel>0) G4cout << "  sigma0=" <<  sigma0<< G4endl;
+
+  const G4double kAnalyticalApproximation= 1.5; 
+ 
+  G4double x = kAnalyticalApproximation/velocity;
+  if (verboseLevel>0) G4cout << "  x=" << x<< G4endl;
+
+
+
+
+  G4double electrIonizationEnergy;                                         
+                                       
+  if ( x<0.035)  
+    {
+      electrIonizationEnergy= 0.75*pi*(std::log(1./(x*x))-1.);  
+    }
+  else 
+    { 
+      if ( x<3.) 
+	{ 
+	  electrIonizationEnergy =std::exp(-2.*x)/(0.031+(0.213*std::pow(x,0.5))+(0.005*x)-(0.069*std::pow(x,3./2.))+(0.324*x*x));
+	}
+     
+      else  
+	{
+	 electrIonizationEnergy =2.*std::exp(-2.*x)/std::pow(x,1.6); }
+    }
+
+  if (verboseLevel>0) G4cout << "  electrIonizationEnergy=" << electrIonizationEnergy<< G4endl;
+
+  G4double hFunction =(electrIonizationEnergy*2.)/(tetaK*std::pow(velocity,3)); //hFunction represents the correction for polarization effet
+  
+  if (verboseLevel>0) G4cout << "  hFunction=" << hFunction<< G4endl;
+    
+  G4double gFunction = (1.+(9.*velocity)+(31.*velocity*velocity)+(98.*std::pow(velocity,3.))+(12.*std::pow(velocity,4.))+(25.*std::pow(velocity,5.))
+			+(4.2*std::pow(velocity,6.))+(0.515*std::pow(velocity,7.)))/std::pow(1.+velocity,9.); //gFunction represents the correction for binding effet
+  if (verboseLevel>0) G4cout << "  gFunction=" << gFunction<< G4endl;
+ 
+  //-----------------------------------------------------------------------------------------------------------------------------
+
+  G4double sigmaPSS = 1.+(((2.*zIncident)/(screenedzTarget*tetaK))*(gFunction-hFunction)); //describes the perturbed stationnairy state of the affected atomic electon
+  if (verboseLevel>0) G4cout << "  sigmaPSS=" << sigmaPSS<< G4endl;
+  if (verboseLevel>0) G4cout << "  sigmaPSS*tetaK=" << sigmaPSS*tetaK<< G4endl;
+
+  //----------------------------------------------------------------------------------------------------------------------------
+  
+  const G4double cNaturalUnit= 1/fine_structure_const;  // it's the speed of light according to Atomic-Unit-System
+  if (verboseLevel>0) G4cout << "  cNaturalUnit=" << cNaturalUnit<< G4endl;
+  
+  G4double ykFormula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocity/sigmaPSS);
+  if (verboseLevel>0) G4cout << "  ykFormula=" << ykFormula<< G4endl;
+ 
+  G4double relativityCorrection = std::pow((1.+(1.1*ykFormula*ykFormula)),0.5)+ykFormula;// the relativistic correction parameter
+  if (verboseLevel>0) G4cout << "  relativityCorrection=" << relativityCorrection<< G4endl;
+
+  G4double reducedVelocity = velocity*std::pow(relativityCorrection,0.5);  // presents the reduced collision velocity parameter 
+  if (verboseLevel>0) G4cout << "  reducedVelocity=" << reducedVelocity<< G4endl;
+
+  G4double etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)
+                           /(sigmaPSS*tetaK)/(sigmaPSS*tetaK);
+  if (verboseLevel>0) G4cout << "  etaOverTheta2=" << etaOverTheta2<< G4endl;
+
+  // low velocity formula
+  
+  G4double universalFunction = (std::pow(2.,9.)/45.)*std::pow(reducedVelocity/sigmaPSS,8.)*std::pow((1.+(1.72*(reducedVelocity/sigmaPSS)*(reducedVelocity/sigmaPSS))),-4.);// is the reduced universal cross section
+  if (verboseLevel>0) G4cout << "  universalFunction by Brandt 1981 =" << universalFunction<< G4endl;
+
+  // Alternative formula by Rice 1977, closer to tabulated data than the above function from Brandt 1981
+     
+     G4double x_ = etaOverTheta2;
+     if (verboseLevel>0) G4cout << "  x_=" << x_ << G4endl;
+
+     G4double b0 = pow(2.,17)/45.;
+     if (verboseLevel>0) G4cout << "  b0=" << b0 << G4endl;
+
+     G4double b1 = 19*(sigmaPSS*tetaK) - 480./11.;
+     if (verboseLevel>0) G4cout << "  b1=" << b1 << G4endl;
+
+     G4double b2 = (720./7.)*(23*(sigmaPSS*tetaK)*(sigmaPSS*tetaK)/11 - 97*(sigmaPSS*tetaK)/9. + 160./13);
+     if (verboseLevel>0) G4cout << "  b2=" << b2 << G4endl;
+
+     universalFunction = b0*x_*x_*x_*x_*(1+b1*x_+b2*x_*x_);
+
+     if (verboseLevel>0) G4cout << "  universalFunction by Rice 1977 =" << universalFunction<< G4endl;
+
+  
+  if ( etaOverTheta2 < 0.01 )
+  {
+     if (verboseLevel>0) G4cout << "  Notice : FK is computed from low velocity formula" << G4endl;
+     
+  }
+
+  else
+  
+  {
+    
+    if ( etaOverTheta2 > 95 && (sigmaPSS*tetaK) > 0.4 && (sigmaPSS*tetaK) < 1.7 )
+    {
+      // From Rice 1977
+      
+      if (verboseLevel>0) G4cout << "  Notice : FK is computed from high velocity formula" << G4endl;
+
+      if (verboseLevel>0) G4cout << "  sigmaPSS*tetaK=" << sigmaPSS*tetaK << G4endl;
+    
+      G4double C1= tableC1->FindValue(sigmaPSS*tetaK);
+      G4double C2= tableC2->FindValue(sigmaPSS*tetaK);
+      G4double C3= tableC3->FindValue(sigmaPSS*tetaK);
+      if (verboseLevel>0) G4cout << "  C1=" << C1 << G4endl;
+      if (verboseLevel>0) G4cout << "  C2=" << C2 << G4endl;
+      if (verboseLevel>0) G4cout << "  C3=" << C3 << G4endl;
+   
+      G4double etaK = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget);
+      if (verboseLevel>0) G4cout << "  etaK=" << etaK << G4endl;
+
+      G4double etaT = (sigmaPSS*tetaK)*(sigmaPSS*tetaK)*(95.); // at any theta, the largest tabulated etaOverTheta2 is 95
+      if (verboseLevel>0) G4cout << "  etaT=" << etaT << G4endl;
+    
+      G4double fKT = FunctionFK((sigmaPSS*tetaK),95.)*(etaT/(sigmaPSS*tetaK));
+      if (FunctionFK((sigmaPSS*tetaK),95.)<=0.) G4cout << 
+      "*** WARNING : G4ecpssrCrossSection::CalculateCrossSection is unable to interpolate FK function in high velocity region ! ***" << G4endl;
+      if (verboseLevel>0) G4cout << "  FunctionFK=" << FunctionFK((sigmaPSS*tetaK),95.) << G4endl;
+      if (verboseLevel>0) G4cout << "  fKT=" << fKT << G4endl;
+    
+      G4double GK = C2/(4*etaK) + C3/(32*etaK*etaK);
+      if (verboseLevel>0) G4cout << "  GK=" << GK << G4endl;
+      G4double GT = C2/(4*etaT) + C3/(32*etaT*etaT);
+      if (verboseLevel>0) G4cout << "  GT=" << GT << G4endl;
+    
+      G4double DT = fKT - C1*std::log(etaT) + GT;
+      if (verboseLevel>0) G4cout << "  DT=" << DT << G4endl;
+
+      G4double fKK = C1*std::log(etaK) + DT - GK;
+      if (verboseLevel>0) G4cout << "  fKK=" << fKK << G4endl;
+    
+      G4double universalFunction3= fKK/(etaK/tetaK);
+      if (verboseLevel>0) G4cout << "  universalFunction3=" << universalFunction3 << G4endl;
+
+      universalFunction=universalFunction3;
+    
+    }
+    else
+    {
+      // From Rice 1977
+      
+      if (verboseLevel>0) G4cout << "  Notice : FK is computed from INTERPOLATED data" << G4endl;
+     
+      G4double universalFunction2 = FunctionFK((sigmaPSS*tetaK),etaOverTheta2);
+      if (universalFunction2<=0) G4cout << 
+      "*** WARNING : G4ecpssrCrossSection::CalculateCrossSection is unable to interpolate FK function in medium velocity region ! ***" << G4endl;
+      
+      if (verboseLevel>0) G4cout << "  universalFunction2=" << universalFunction2 << " for theta=" << sigmaPSS*tetaK << " and etaOverTheta2=" << etaOverTheta2 << G4endl;
+      
+      universalFunction=universalFunction2;
+    }
+        
+  }
+  
+  //----------------------------------------------------------------------------------------------------------------------
+
+  G4double sigmaPSSR = (sigma0/(sigmaPSS*tetaK))*universalFunction; //sigmaPSSR is the straight-line K-shell ionization cross section
+  if (verboseLevel>0) G4cout << "  sigmaPSSR=" << sigmaPSSR<< G4endl;
+  
+  //-----------------------------------------------------------------------------------------------------------------------
+
+  G4double pssDeltaK = (4./(systemMass*sigmaPSS*tetaK))*(sigmaPSS/velocity)*(sigmaPSS/velocity);
+  if (verboseLevel>0) G4cout << "  pssDeltaK=" << pssDeltaK<< G4endl;
+
+  G4double energyLoss = std::pow(1-pssDeltaK,0.5); //energyLoss incorporates the straight-line energy-loss 
+  if (verboseLevel>0) G4cout << "  energyLoss=" << energyLoss<< G4endl;
+
+  G4double energyLossFunction = (std::pow(2.,-9)/8.)*((((9.*energyLoss)-1.)*std::pow(1.+energyLoss,9.))+(((9.*energyLoss)+1.)*std::pow(1.-energyLoss,9.)));//energy loss function 
+  if (verboseLevel>0) G4cout << "  energyLossFunction=" <<  energyLossFunction<< G4endl;
+
+  //----------------------------------------------------------------------------------------------------------------------------------------------
+
+  G4double coulombDeflection = (4.*pi*zIncident/systemMass)*std::pow(tetaK*sigmaPSS,-2.)*std::pow(velocity/sigmaPSS,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter 
+ 
+  G4double cParameter = 2.*coulombDeflection/(energyLoss*(energyLoss+1.));
+  
+  if (verboseLevel>0) G4cout << "  cParameter=" << cParameter<< G4endl;
+
+  G4double coulombDeflectionFunction = 9.*ExpIntFunction(10,cParameter);                         //this function describes Coulomb-deflection effect
+  if (verboseLevel>0) G4cout << "  ExpIntFunction(10,cParameter) =" << ExpIntFunction(10,cParameter) << G4endl;
+  if (verboseLevel>0) G4cout << "  coulombDeflectionFunction =" << coulombDeflectionFunction << G4endl;
+
+  //--------------------------------------------------------------------------------------------------------------------------------------------------
+ 
+  G4double crossSection =  0;
+  
+  crossSection = energyLossFunction* coulombDeflectionFunction*sigmaPSSR;  //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS)
+                                                                           //and it's reduced by the energy-loss(E),the Coulomb deflection(C),
+                                                                           //and the relativity(R) effects
+
+  //--------------------------------------------------------------------------------------------------------------------------------------------------   
+
+  if (crossSection >= 0) {
+    return crossSection;
+  }
+  else {return 0;}
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ecpssrKCrossSection::FunctionFK(G4double k, G4double theta) 							  
+{
+
+  G4double sigma = 0.;
+  G4double valueT1 = 0;
+  G4double valueT2 = 0;
+  G4double valueE21 = 0;
+  G4double valueE22 = 0;
+  G4double valueE12 = 0;
+  G4double valueE11 = 0;
+  G4double xs11 = 0;   
+  G4double xs12 = 0; 
+  G4double xs21 = 0; 
+  G4double xs22 = 0; 
+
+  // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values 
+  // (in particular for FK computation at 95 for high velocity formula)
+  
+  if (
+  theta==9.5e-2 ||
+  theta==9.5e-1 ||
+  theta==9.5e+00 ||
+  theta==9.5e+01
+  ) theta=theta-1e-12;
+
+  if (
+  theta==1.e-2 ||
+  theta==1.e-1 ||
+  theta==1.e+00 ||
+  theta==1.e+01
+  ) theta=theta+1e-12;
+
+  // END PROTECTION
+  
+  {
+    std::vector<double>::iterator t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
+    std::vector<double>::iterator t1 = t2-1;
+ 
+    std::vector<double>::iterator e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
+    std::vector<double>::iterator e11 = e12-1; 
+	  
+    std::vector<double>::iterator e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
+    std::vector<double>::iterator e21 = e22-1;
+	  	
+    valueT1  =*t1;
+    valueT2  =*t2;
+    valueE21 =*e21;
+    valueE22 =*e22;
+    valueE12 =*e12;
+    valueE11 =*e11;
+
+    xs11 = FKData[valueT1][valueE11];
+    xs12 = FKData[valueT1][valueE12];
+    xs21 = FKData[valueT2][valueE21];
+    xs22 = FKData[valueT2][valueE22];
+
+/*
+verboseLevel=1;
+
+if (verboseLevel>0) G4cout << "x1= " << valueT1 << G4endl;
+if (verboseLevel>0) G4cout << " vector of y for x1" << G4endl;
+for_each (aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), print);
+
+if (verboseLevel>0) G4cout << G4endl;
+
+if (verboseLevel>0) G4cout << "x2= " << valueT2 << G4endl;
+if (verboseLevel>0) G4cout << " vector of y for x2" << G4endl;
+for_each (aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), print);
+
+if (verboseLevel>0) G4cout << G4endl;
+
+if (verboseLevel>0) G4cout 
+<< "  " 
+<< valueT1 << " "
+<< valueT2 << " "
+<< valueE11 << " "
+<< valueE12 << " "
+<< valueE21<< " "
+<< valueE22 << " " 
+<< xs11 << " " 
+<< xs12 << " " 
+<< xs21 << " " 
+<< xs22 << " " 
+<< G4endl;
+//verboseLevel=0;
+*/
+
+}
+     
+  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
+  
+  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
+     
+  if (xsProduct != 0.)
+  {
+    sigma = QuadInterpolator(  valueE11, valueE12, 
+    			       valueE21, valueE22, 
+			       xs11, xs12, 
+			       xs21, xs22, 
+			       valueT1, valueT2, 
+			       k, theta );
+  }
+
+  return sigma;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ecpssrKCrossSection::LinLogInterpolate(G4double e1, 
+						        G4double e2, 
+						        G4double e, 
+						        G4double xs1, 
+						        G4double xs2)
+{
+  G4double d1 = std::log(xs1);
+  G4double d2 = std::log(xs2);
+  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
+  return value;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ecpssrKCrossSection::LogLogInterpolate(G4double e1, 
+						        G4double e2, 
+						        G4double e, 
+						        G4double xs1, 
+						        G4double xs2)
+{
+  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
+  G4double b = std::log10(xs2) - a*std::log10(e2);
+  G4double sigma = a*std::log10(e) + b;
+  G4double value = (std::pow(10.,sigma));
+  return value;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ecpssrKCrossSection::QuadInterpolator(G4double e11, G4double e12, 
+						       G4double e21, G4double e22, 
+						       G4double xs11, G4double xs12, 
+						       G4double xs21, G4double xs22, 
+						       G4double t1, G4double t2, 
+						       G4double t, G4double e)
+{
+// Log-Log
+/*
+  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
+  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
+  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
+*/
+
+// Lin-Log
+  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
+  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
+  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
+  return value;
+}
+
+
+
+
+
Index: trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrLiCrossSection.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrLiCrossSection.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4ecpssrLiCrossSection.cc	(revision 1197)
@@ -0,0 +1,982 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//$Id: G4ecpssrLiCrossSection.cc,v 1.7 2009/11/11 09:14:53 mantero Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+// Author: Haifa Ben Abdelouahed
+//         
+//
+// History:
+// -----------
+//  23 Apr 2008   H. Ben Abdelouahed   1st implementation
+//  28 Apr 2008   MGP        Major revision according to a design iteration
+//  29 Apr 2009   ALF Updated Desing for Integration
+//  02 May 2009   ALF + Haifa L1,L2,L3 Extensions
+//  11 Nov 2009   ALF code cleaning for the Dec release
+//
+// -------------------------------------------------------------------
+// Class description:
+// Low Energy Electromagnetic Physics, Cross section, p and alpha ionisation, L shell
+// Further documentation available from http://www.ge.infn.it/geant4/lowE
+// -------------------------------------------------------------------
+
+
+#include "globals.hh"
+#include "G4ecpssrLiCrossSection.hh"
+#include "G4AtomicTransitionManager.hh"
+#include "G4NistManager.hh"
+#include "G4Proton.hh"
+#include "G4Alpha.hh"
+#include <math.h>
+
+G4ecpssrLiCrossSection::G4ecpssrLiCrossSection()
+{ 
+
+
+ // Storing FLi data needed for 0.2 to 3.0  velocities region
+
+    char *path = getenv("G4LEDATA");
+ 
+    if (!path)
+    G4Exception("G4ecpssrLCrossSection::CalculateCrossSection: G4LEDDATA environment variable not set");
+
+     std::ostringstream fileName1;
+     std::ostringstream fileName2;
+    fileName1 << path << "/pixe/uf/FL1.dat";
+    fileName2 << path << "/pixe/uf/FL2.dat";
+   
+    std::ifstream FL1(fileName1.str().c_str());
+    std::ifstream FL2(fileName1.str().c_str());
+   
+   
+    if (!FL1) G4Exception("G4ecpssrLCrossSection::CalculateCrossSection: error opening FL1 data file");
+    if (!FL2) G4Exception("G4ecpssrLCrossSection::CalculateCrossSection: error opening FL2 data file");
+  
+ dummyVec.push_back(0.);
+
+    while(!FL1.eof())
+    {
+	double x1;
+	double y1;
+	
+	FL1>>x1>>y1;
+
+	//  Mandatory vector initialization
+        if (x1 != dummyVec.back()) 
+        { 
+          dummyVec.push_back(x1); 
+          aVecMap[x1].push_back(-1.);
+        }
+	  
+        FL1>>FL1Data[x1][y1];
+
+        if (y1 != aVecMap[x1].back()) aVecMap[x1].push_back(y1);
+    }
+ while(!FL2.eof())
+    {
+	double x2;
+	double y2;
+	
+	FL2>>x2>>y2;
+
+	//  Mandatory vector initialization
+        if (x2 != dummyVec.back()) 
+        { 
+          dummyVec.push_back(x2); 
+          aVecMap[x2].push_back(-1.);
+        }
+	  
+        FL2>>FL2Data[x2][y2];
+
+        if (y2 != aVecMap[x2].back()) aVecMap[x2].push_back(y2);    
+    }
+
+
+}
+
+G4ecpssrLiCrossSection::~G4ecpssrLiCrossSection()
+{ }
+
+//---------------------------------this "ExpIntFunction" function allows fast evaluation of the n order exponential integral function En(x)------
+
+G4double G4ecpssrLiCrossSection::ExpIntFunction(G4int n,G4double x)
+
+{
+  G4int i;
+  G4int ii;
+  G4int nm1;
+  G4double a;
+  G4double b;
+  G4double c;
+  G4double d;
+  G4double del;
+  G4double fact;
+  G4double h;
+  G4double psi;
+  G4double ans = 0;
+  const G4double euler= 0.5772156649;
+  const G4int maxit= 100;
+  const G4double fpmin = 1.0e-30;
+  const G4double eps = 1.0e-7;
+  nm1=n-1;
+  if (n<0 || x<0.0 || (x==0.0 && (n==0 || n==1)))
+  G4cout << "bad arguments in ExpIntFunction" << G4endl;
+  else {
+       if (n==0) ans=exp(-x)/x;
+        else {
+           if (x==0.0) ans=1.0/nm1;
+              else {
+                   if (x > 1.0) {
+                         	b=x+n;
+                        	c=1.0/fpmin;
+                        	d=1.0/b;
+				h=d;
+				for (i=1;i<=maxit;i++) {
+				  a=-i*(nm1+i);
+				  b +=2.0;
+				  d=1.0/(a*d+b);
+				  c=b+a/c;
+				  del=c*d;
+				  h *=del;
+				      if (fabs(del-1.0) < eps) {
+					ans=h*exp(-x);
+					return ans;
+				      }
+				}
+		   } else {
+		     ans = (nm1!=0 ? 1.0/nm1 : -log(x)-euler);
+		     fact=1.0;
+		     for (i=1;i<=maxit;i++) {
+		       fact *=-x/i;
+		       if (i !=nm1) del = -fact/(i-nm1);
+		       else {
+			 psi = -euler;
+			 for (ii=1;ii<=nm1;ii++) psi +=1.0/ii;
+			 del=fact*(-log(x)+psi);
+		       }
+		       ans += del;
+		       if (fabs(del) < fabs(ans)*eps) return ans;
+		     }
+		   }
+	      }
+	}
+  }
+return ans;
+}
+//-----------------------------------------------------------------------------------------------------------
+
+ 
+G4double G4ecpssrLiCrossSection::CalculateL1CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
+ 
+ //this L-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//					                
+
+{
+
+  G4NistManager* massManager = G4NistManager::Instance();   
+
+  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
+
+
+  G4int zIncident = 0; 
+  G4Proton* aProtone = G4Proton::Proton();
+  G4Alpha* aAlpha = G4Alpha::Alpha();
+
+  if (massIncident == aProtone->GetPDGMass() )
+  {
+
+   zIncident = (G4int)((aProtone->GetPDGCharge())/eplus); 
+
+   //G4cout << "zincident:" << zIncident << G4endl;
+  }
+  else
+    {
+      if (massIncident == aAlpha->GetPDGMass())
+	{
+	  
+	  zIncident  =(G4int) ((aAlpha->GetPDGCharge())/eplus); 
+	  
+	  //G4cout << "zincident:" << zIncident << G4endl;
+	}
+      else
+	{ 
+	  G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
+	  massIncident =0.;
+	}
+    }
+
+
+  G4double l1BindingEnergy = transitionManager->Shell(zTarget,1)->BindingEnergy();
+  
+
+
+
+
+         
+  G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
+ 
+  G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//the mass of the system (projectile, target) 
+
+  const G4double zlshell= 4.15;
+
+  G4double screenedzTarget = zTarget-zlshell;         // screenedzTarget is the screened nuclear charge of the target
+
+  const G4double rydbergMeV= 13.6056923e-6;
+
+  const G4double nl= 2.;        // nl is the quantum number of the L shell 
+     
+  G4double tetal1 = (l1BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV);  //tetal1 denotes the reduced L1-shell-binding-energy of the electron
+
+
+
+
+
+
+  const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;        
+  
+  G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*pow(screenedzTarget,-4.);     //sigma0 is the initial cross section of L shell at stable state
+
+  //---------------------------------------------------------------------------------------------------------------------
+
+  G4double velocityl1 = CalculateVelocity(1, zTarget, massIncident, energyIncident);   //is the scaled velocity parameter of the system  
+
+
+
+  
+  //---------------------------------------------------------------------------------------------------------------------
+  
+  const G4double l1AnalyticalApproximation= 1.5; 
+
+
+
+  G4double x1 = nl*l1AnalyticalApproximation/velocityl1;
+
+
+
+
+
+  //-----------------------------------------x of l1 sub shell--------------------------------------
+ 
+ G4double electrIonizationEnergyl1;
+                                       
+                                       
+  if ( x1<0.035 && x1>= 0.)  
+    {
+      electrIonizationEnergyl1= 0.75*pi*(log(1./(x1*x1))-1.);  
+    }
+  else 
+    { 
+      if ( x1<3.&& x1>=0.035) 
+	{ 
+	  electrIonizationEnergyl1 =exp(-2.*x1)/(0.031+(0.213*pow(x1,0.5))+(0.005*x1)-(0.069*pow(x1,3./2.))+(0.324*x1*x1));
+	}
+     
+      else  
+	{
+	  if ( x1<=11.&& x1>=3.) {
+	    electrIonizationEnergyl1 =2.*exp(-2.*x1)/pow(x1,1.6);
+	  }
+	  else {
+	    electrIonizationEnergyl1 =0.;
+	      }
+	}
+    }
+
+
+ 
+ 
+
+  //-------------------------------------------------------- h and g functions for l1 -------------------------------------------------
+
+
+  G4double hFunctionl1 =(electrIonizationEnergyl1*2.*nl)/(tetal1*pow(velocityl1,3)); //hFunction represents the correction for polarization effet
+    
+ 
+  G4double gFunctionl1 = (1.+(9.*velocityl1)+(31.*velocityl1*velocityl1)+(49.*pow(velocityl1,3.))+(162.*pow(velocityl1,4.))+(63.*pow(velocityl1,5.))
+			  +(18.*pow(velocityl1,6.))+(1.97*pow(velocityl1,7.)))/pow(1.+velocityl1,9.); //gFunction represents the correction for binding effet
+
+
+
+
+  //-----------------------------------------------------------------------------------------------------------------------------
+
+  G4double sigmaPSS_l1 = 1.+(((2.*zIncident)/(screenedzTarget*tetal1))*(gFunctionl1-hFunctionl1)); //describes the perturbed stationnairy state of the affected atomic electon
+ 
+ 
+
+ 
+
+
+
+ //----------------------------------------------------------------------------------------------------------------------------
+  //const G4double cNaturalUnit= 1/fine_structure_const;  // it's the speed of light according to Atomic-Unit-System
+  //--------------------------------------------------------------------------------------------------------------
+
+  //G4double yl1Formula=0.4*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(nl*(velocityl1/sigmaPSS_l1));
+
+ 
+
+  //-----------------------------------------------Relativity effect correction L1 -------------------------------------------------------------
+
+  //G4double relativityCorrectionl1 = pow((1.+(1.1*yl1Formula*yl1Formula)),0.5)+yl1Formula;// the relativistic correction parameter
+  //G4double reducedVelocityl1 = velocityl1*pow(relativityCorrectionl1,0.5);  // presents the reduced collision velocity parameter 
+
+
+
+
+
+
+
+  //-------------------------------------------------------------------------------------------------------------------
+   //------------------------------------------------------------UNIVERSAL FUNCTION ---------------------------------
+  //------------------------------------------------------------------------------------------------------------
+                                                  // is the reduced universal cross section
+
+  G4double universalFunction_l1 ;
+
+ 
+ 
+    //-------------------------------------------------------------------------------------------------------------------
+   //------------------------------------------------------------ LIMITS OF ECPSSR MODEL ---------------------------
+  //------------------------------------------------------------------------------------------------------------
+
+
+  G4double L1etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)/(sigmaPSS_l1*tetal1)/(sigmaPSS_l1*tetal1);
+
+
+     universalFunction_l1 = FunctionFL1((sigmaPSS_l1*tetal1), L1etaOverTheta2);
+        
+
+
+
+     //-----------------------------------------------------------------------------------------------------------------------
+    //---------------------------------------------------------------------------   PSSR Li CROSS SECTION  -------------
+   //-----------------------------------------------------------------------------------------------------------------------
+
+   G4double sigmaPSSR_l1 = (sigma0/(sigmaPSS_l1*tetal1))*universalFunction_l1; //sigmaPSSR is the straight-line L-shell ionization cross section
+
+
+
+   //----------------------------------------------------------------------------------------------------------------------
+  
+  G4double pssDeltal1 = (4./(systemMass*sigmaPSS_l1*tetal1))*(sigmaPSS_l1/velocityl1)*(sigmaPSS_l1/velocityl1);
+
+
+
+  G4double energyLossl1 = pow(1-pssDeltal1,0.5); //energyLoss incorporates the straight-line energy-loss 
+
+
+
+
+
+     //-----------------------------------------------------------------------------------------------------------------------
+    //----------------------------------------------------------------------------ENERGY LOSS CORRECTION-------------
+   //-----------------------------------------------------------------------------------------------------------------------
+
+/* 
+ G4double energyLossFunction_L1 = (pow(2.,-9)/8.)*((((9.*energyLossl1)-1.)*pow(1.+energyLossl1,9.))+(((9.*energyLossl1)+1.)*pow(1.-energyLossl1,9.))); 
+  G4double energyLossFunction_L2 = (pow(2.,-11)/10.)*((((11.*energyLossl2)-1.)*pow(1.+energyLossl2,11.))+(((11.*energyLossl2)+1.)*pow(1.-energyLossl2,11.)));
+  G4double energyLossFunction_L3 = (pow(2.,-11)/10.)*((((11.*energyLossl3)-1.)*pow(1.+energyLossl3,11.))+(((11.*energyLossl3)+1.)*pow(1.-energyLossl3,11.)));
+*/
+
+
+     //-----------------------------------------------------------------------------------------------------------------------
+    //----------------------------------------------------------------------------COULOMB DEFLECTION CORRECTION-------------
+   //-----------------------------------------------------------------------------------------------------------------------
+
+
+  G4double coulombDeflectionl1 = (4.*pi*zIncident/systemMass)*pow(tetal1*sigmaPSS_l1,-2.)*pow(velocityl1/sigmaPSS_l1,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter 
+
+
+
+
+  G4double cParameterl1 = 2.*coulombDeflectionl1/(energyLossl1*(energyLossl1+1.));
+
+
+
+  G4double coulombDeflectionFunction_l1 = 9.*ExpIntFunction(10,cParameterl1);                         //this function describes Coulomb-deflection effect
+
+
+
+      //--------------------------------------------------------------------------------------------------------------------------------------------------
+     //-----------------------------------------------------------------------------------------------------------------------
+    //---------------------------------------------------------------------------   ECPSSR Li CROSS SECTION  -------------
+   //-----------------------------------------------------------------------------------------------------------------------
+  //--------------------------------------------------------------------------------------------------------------------------------------------------
+ 
+  /*
+   G4double crossSection_L1 = energyLossFunction_L1 * coulombDeflectionFunction_l1 * sigmaPSSR_l1;  //this ECPSSR cross section is estimated at perturbed-stationnairy-state(PSS)
+   G4double crossSection_L2 = energyLossFunction_L2 * coulombDeflectionFunction_l2 * sigmaPSSR_l2;         //and it's reduced by the energy-loss(E),the Coulomb deflection(C),
+   G4double crossSection_L3 = energyLossFunction_L3 * coulombDeflectionFunction_l3 * sigmaPSSR_l3;              //and the relativity(R) effects
+  */
+
+   G4double crossSection_L1 =  coulombDeflectionFunction_l1 * sigmaPSSR_l1; 
+
+ 
+ if (crossSection_L1 >= 0) {
+    return crossSection_L1;
+  }
+  else {return 0;}
+}
+
+G4double G4ecpssrLiCrossSection::CalculateL2CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
+ 
+ //this L-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//					                
+
+{
+
+
+  G4NistManager* massManager = G4NistManager::Instance();   
+
+  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
+
+
+  G4int zIncident = 0; 
+  G4Proton* aProtone = G4Proton::Proton();
+  G4Alpha* aAlpha = G4Alpha::Alpha();
+
+  if (massIncident == aProtone->GetPDGMass() )
+  {
+
+   zIncident =(G4int) ((aProtone->GetPDGCharge())/eplus); 
+
+   //   G4cout << "zincident:" << zIncident << G4endl;
+    }
+  else
+    {
+      if (massIncident == aAlpha->GetPDGMass())
+	{
+	  
+	  zIncident  = (G4int) ((aAlpha->GetPDGCharge())/eplus); 
+	  
+	  //	  G4cout << "zincident:" << zIncident << G4endl;
+	}
+      else
+	{ 
+	  G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
+	  massIncident =0.;
+	}
+    }
+
+  G4double l2BindingEnergy = transitionManager->Shell(zTarget,2)->BindingEnergy();
+
+  G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
+  G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//the mass of the system (projectile, target) 
+  const G4double zlshell= 4.15;
+  G4double screenedzTarget = zTarget-zlshell;         // screenedzTarget is the screened nuclear charge of the target
+  const G4double rydbergMeV= 13.6056923e-6;
+  const G4double nl= 2.;        // nl is the quantum number of the L shell 
+
+  G4double tetal2 = (l2BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); 
+
+  const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;        
+  G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*pow(screenedzTarget,-4.);     //sigma0 is the initial cross section of L shell at stable state
+
+  G4double velocityl2 = CalculateVelocity(2,  zTarget, massIncident, energyIncident);
+
+  const G4double l2AnalyticalApproximation= 1.25;
+
+  G4double x2 = nl*l2AnalyticalApproximation/velocityl2;
+
+  //---------------------------------------------------------------x of l2 sub shell------------------------------------------
+ 
+G4double electrIonizationEnergyl2;
+                                       
+                                       
+  if ( x2<0.035 && x2>= 0.)  
+    {
+      electrIonizationEnergyl2= 0.75*pi*(log(1./(x2*x2))-1.);  
+    }
+  else 
+    { 
+      if ( x2<3. && x2 >=0.035) 
+	{ 
+	  electrIonizationEnergyl2 =exp(-2.*x2)/(0.031+(0.213*pow(x2,0.5))+(0.005*x2)-(0.069*pow(x2,3./2.))+(0.324*x2*x2));
+	}
+     
+      else  
+	{
+	  
+	  if ( x2<=11.&& x2>=3.) {
+	    electrIonizationEnergyl2 =2.*exp(-2.*x2)/pow(x2,1.6);
+	  }
+	  else {
+	    electrIonizationEnergyl2=0.;
+	  }
+	  
+	}
+    }
+  
+  //----------------------------------------------------------------------------- h and g functions for l2  ----------------------------
+  
+  
+  G4double hFunctionl2 =(electrIonizationEnergyl2*2.*nl)/(tetal2*pow(velocityl2,3)); //hFunction represents the correction for polarization effet
+  
+  G4double gFunctionl2 = (1.+(10.*velocityl2)+(45.*velocityl2*velocityl2)+(102.*pow(velocityl2,3.))+(331.*pow(velocityl2,4.))+(6.7*pow(velocityl2,5.))
+			  +(58.*pow(velocityl2,6.))+(7.8*pow(velocityl2,7.))+ (0.888*pow(velocityl2,8.)) )/pow(1.+velocityl2,10.); //gFunction represents the correction for binding effet
+  
+  
+  G4double sigmaPSS_l2 = 1.+(((2.*zIncident)/(screenedzTarget*tetal2))*(gFunctionl2-hFunctionl2)); 
+  
+  //----------------------------------------------------------------------------------------------------------------------------
+  //const G4double cNaturalUnit= 1/fine_structure_const;  // it's the speed of light according to Atomic-Unit-System
+  //--------------------------------------------------------------------------------------------------------------
+  
+  //G4double yl2Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl2/sigmaPSS_l2);
+  
+  //-----------------------------------------------Relativity effect correction L2 --------------------------------------------------------------
+  
+  //G4double relativityCorrectionl2 = pow((1.+(1.1*yl2Formula*yl2Formula)),0.5)+yl2Formula;// the relativistic correction parameter
+  //  G4double reducedVelocityl2 = velocityl2*pow(relativityCorrectionl2,0.5);  // presents the reduced collision velocity parameter 
+  
+  G4double L2etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)/(sigmaPSS_l2*tetal2)/(sigmaPSS_l2*tetal2);
+  
+  
+  G4double universalFunction_l2 ;
+  
+  
+  universalFunction_l2 = FunctionFL2((sigmaPSS_l2*tetal2), L2etaOverTheta2);
+  
+  G4double sigmaPSSR_l2 = (sigma0/(sigmaPSS_l2*tetal2))*universalFunction_l2;
+  G4double pssDeltal2 = (4./(systemMass*sigmaPSS_l2*tetal2))*(sigmaPSS_l2/velocityl2)*(sigmaPSS_l2/velocityl2);
+  G4double energyLossl2 = pow(1-pssDeltal2,0.5);  
+  
+  G4double coulombDeflectionl2 = (4.*pi*zIncident/systemMass)*pow(tetal2*sigmaPSS_l2,-2.)*pow(velocityl2/sigmaPSS_l2,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter 
+  G4double cParameterl2 = 2.*coulombDeflectionl2/(energyLossl2*(energyLossl2+1.));
+  G4double coulombDeflectionFunction_l2 = 11.*ExpIntFunction(12,cParameterl2);
+  
+  G4double crossSection_L2 = coulombDeflectionFunction_l2 * sigmaPSSR_l2 ;         
+  
+  
+  if (crossSection_L2 >= 0) {
+    return crossSection_L2;
+  }
+  else {return 0;}
+  
+  
+}
+
+
+G4double G4ecpssrLiCrossSection::CalculateL3CrossSection(G4int zTarget,G4double massIncident, G4double energyIncident)
+ 
+ //this L-CrossSection calculation method is done according to W.Brandt and G.Lapicki, Phys.Rev.A23(1981)//					                
+
+{
+
+
+  G4NistManager* massManager = G4NistManager::Instance();   
+
+  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
+
+
+  G4int zIncident = 0; 
+  G4Proton* aProtone = G4Proton::Proton();
+  G4Alpha* aAlpha = G4Alpha::Alpha();
+
+  if (massIncident == aProtone->GetPDGMass() )
+  {
+
+   zIncident = (G4int) ((aProtone->GetPDGCharge())/eplus); 
+
+   //   G4cout << "zincident:" << zIncident << G4endl;
+    }
+  else
+    {
+      if (massIncident == aAlpha->GetPDGMass())
+	{
+	  
+	  zIncident  =(G4int) ((aAlpha->GetPDGCharge())/eplus); 
+	  
+	  //	  G4cout << "zincident:" << zIncident << G4endl;
+	}
+      else
+	{ 
+	  G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
+	  massIncident =0.;
+	}
+    }
+
+  G4double l3BindingEnergy = transitionManager->Shell(zTarget,3)->BindingEnergy();
+        
+  G4double massTarget = (massManager->GetAtomicMassAmu(zTarget))*amu_c2;
+ 
+  G4double systemMass =((massIncident*massTarget)/(massIncident+massTarget))/electron_mass_c2;//the mass of the system (projectile, target) 
+
+  const G4double zlshell= 4.15;
+
+  G4double screenedzTarget = zTarget-zlshell;         // screenedzTarget is the screened nuclear charge of the target
+
+  const G4double rydbergMeV= 13.6056923e-6;
+
+  const G4double nl= 2.;        // nl is the quantum number of the L shell 
+     
+  G4double tetal3 = (l3BindingEnergy*nl*nl)/((screenedzTarget*screenedzTarget)*rydbergMeV); 
+
+
+
+  const G4double bohrPow2Barn=(Bohr_radius*Bohr_radius)/barn ;        
+  
+  G4double sigma0 = 8.*pi*(zIncident*zIncident)*bohrPow2Barn*pow(screenedzTarget,-4.);     //sigma0 is the initial cross section of L shell at stable state
+
+  //---------------------------------------------------------------------------------------------------------------------
+
+  G4double velocityl3 = CalculateVelocity(3, zTarget, massIncident, energyIncident);
+
+  const G4double l3AnalyticalApproximation= 1.25;
+
+  G4double x3 = nl*l3AnalyticalApproximation/velocityl3;
+
+
+  //--------------------------------------------------------------------- x of l3 sub shell---------------------------------
+
+
+ G4double electrIonizationEnergyl3;
+                                     
+                                       
+  if ( x3<0.035 && x3>=0. )  
+    {
+      electrIonizationEnergyl3= 0.75*pi*(log(1./(x3*x3))-1.);  
+    }
+  else 
+    { 
+      if ( x3<3. && x3 >= 0.035) 
+	{ 
+	  electrIonizationEnergyl3 =exp(-2.*x3)/(0.031+(0.213*pow(x3,0.5))+(0.005*x3)-(0.069*pow(x3,3./2.))+(0.324*x3*x3));
+	}
+     
+      else  
+	{
+	  if ( x3<=11.&& x3>=3.) {
+	    electrIonizationEnergyl3 =2.*exp(-2.*x3)/pow(x3,1.6);
+	  }
+	  else {
+            electrIonizationEnergyl3=0.;
+	  }
+	  
+	}
+    }
+
+ //------------------------------------------------------------ h and g function  for l3  ---------------------------------------------
+
+
+  G4double hFunctionl3 =(electrIonizationEnergyl3*2.*nl)/(tetal3*pow(velocityl3,3)); //hFunction represents the correction for polarization effet
+    
+ 
+  G4double gFunctionl3 = (1.+(10.*velocityl3)+(45.*velocityl3*velocityl3)+(102.*pow(velocityl3,3.))+(331.*pow(velocityl3,4.))+(6.7*pow(velocityl3,5.))
+			  +(58.*pow(velocityl3,6.))+(7.8*pow(velocityl3,7.))+ (0.888*pow(velocityl3,8.)) )/pow(1.+velocityl3,10.); //gFunction represents the correction for binding effet
+  //-----------------------------------------------------------------------------------------------------------------------------
+
+  G4double sigmaPSS_l3 = 1.+(((2.*zIncident)/(screenedzTarget*tetal3))*(gFunctionl3-hFunctionl3)); 
+
+
+  //----------------------------------------------------------------------------------------------------------------------------
+  //const G4double cNaturalUnit= 1/fine_structure_const;  // it's the speed of light according to Atomic-Unit-System
+  //--------------------------------------------------------------------------------------------------------------
+
+  //G4double yl3Formula=0.15*(screenedzTarget/cNaturalUnit)*(screenedzTarget/cNaturalUnit)/(velocityl3/sigmaPSS_l3);
+
+ //-----------------------------------------------Relativity effect correction L3 --------------------------------------------------------------
+
+ //G4double relativityCorrectionl3 = pow((1.+(1.1*yl3Formula*yl3Formula)),0.5)+yl3Formula;// the relativistic correction parameter
+ //G4double reducedVelocityl3 = velocityl3*pow(relativityCorrectionl3,0.5);  // presents the reduced collision velocity parameter 
+
+  G4double universalFunction_l3 ;
+
+  G4double L3etaOverTheta2 = (energyIncident*electron_mass_c2)/(massIncident*rydbergMeV*screenedzTarget*screenedzTarget)/(sigmaPSS_l3*tetal3)/(sigmaPSS_l3*tetal3);
+  
+  universalFunction_l3 = 2*FunctionFL2((sigmaPSS_l3*tetal3), L3etaOverTheta2);
+  
+  
+  G4double sigmaPSSR_l3 = (sigma0/(sigmaPSS_l3*tetal3))*universalFunction_l3;
+  
+  G4double pssDeltal3 = (4./(systemMass*sigmaPSS_l3*tetal3))*(sigmaPSS_l3/velocityl3)*(sigmaPSS_l3/velocityl3);
+
+  G4double energyLossl3 = pow(1-pssDeltal3,0.5); 
+
+  G4double coulombDeflectionl3 = (4.*pi*zIncident/systemMass)*pow(tetal3*sigmaPSS_l3,-2.)*pow(velocityl3/sigmaPSS_l3,-3.)*(zTarget/screenedzTarget); //incorporates Coulomb deflection parameter 
+
+  G4double cParameterl3 = 2.*coulombDeflectionl3/(energyLossl3*(energyLossl3+1.));
+  
+  G4double coulombDeflectionFunction_l3 = 11.*ExpIntFunction(12,cParameterl3);
+  
+  G4double crossSection_L3 = coulombDeflectionFunction_l3 * sigmaPSSR_l3;              
+  
+  if (crossSection_L3 >= 0) {
+    return crossSection_L3;
+  }
+  else {return 0;}
+  
+}
+
+
+
+G4double G4ecpssrLiCrossSection::CalculateVelocity(G4int subShell, G4int zTarget, G4double massIncident,  G4double energyIncident) 
+			                     
+{  
+
+  G4AtomicTransitionManager*  transitionManager =  G4AtomicTransitionManager::Instance();
+
+  G4double liBindingEnergy = transitionManager->Shell(zTarget,subShell)->BindingEnergy();
+
+  
+  G4Proton* aProtone = G4Proton::Proton();
+  G4Alpha* aAlpha = G4Alpha::Alpha();
+  
+  if (!(massIncident == aProtone->GetPDGMass() || massIncident == aAlpha->GetPDGMass()))
+    {
+      G4cout << "we can treat only Proton or Alpha incident particles " << G4endl;
+      return 0;
+    }
+  
+  
+  const G4double zlshell= 4.15;      
+  
+  G4double screenedzTarget = zTarget- zlshell;                                  
+  
+  const G4double rydbergMeV= 13.6e-6;
+  
+  const G4double nl= 2.;        // nl is the quantum number of the L shell 
+  
+  G4double tetali = (liBindingEnergy*nl*nl)/(screenedzTarget*screenedzTarget*rydbergMeV);     
+
+   G4double velocity =(2.*nl/(tetali*screenedzTarget))*pow(((energyIncident*electron_mass_c2)/(massIncident*rydbergMeV)),0.5);
+  
+  return velocity;
+}
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+
+G4double G4ecpssrLiCrossSection::FunctionFL1(G4double k, G4double theta) 							  
+{
+
+  G4double sigma = 0.;
+  G4double valueT1 = 0;
+  G4double valueT2 = 0;
+  G4double valueE21 = 0;
+  G4double valueE22 = 0;
+  G4double valueE12 = 0;
+  G4double valueE11 = 0;
+  G4double xs11 = 0;   
+  G4double xs12 = 0; 
+  G4double xs21 = 0; 
+  G4double xs22 = 0; 
+
+  // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values 
+  // (in particular for FK computation at 95 for high velocity formula)
+  
+  if (
+  theta==9.5e-2 ||
+  theta==9.5e-1 ||
+  theta==9.5e+00 ||
+  theta==9.5e+01
+  ) theta=theta-1e-12;
+
+  if (
+  theta==1.e-2 ||
+  theta==1.e-1 ||
+  theta==1.e+00 ||
+  theta==1.e+01
+  ) theta=theta+1e-12;
+
+  // END PROTECTION
+  
+  {
+    std::vector<double>::iterator t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
+    std::vector<double>::iterator t1 = t2-1;
+ 
+    std::vector<double>::iterator e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
+    std::vector<double>::iterator e11 = e12-1; 
+	  
+    std::vector<double>::iterator e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
+    std::vector<double>::iterator e21 = e22-1;
+	  	
+    valueT1  =*t1;
+    valueT2  =*t2;
+    valueE21 =*e21;
+    valueE22 =*e22;
+    valueE12 =*e12;
+    valueE11 =*e11;
+
+    xs11 = FL1Data[valueT1][valueE11];
+    xs12 = FL1Data[valueT1][valueE12];
+    xs21 = FL1Data[valueT2][valueE21];
+    xs22 = FL1Data[valueT2][valueE22];
+
+}
+     
+  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
+  
+  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
+     
+  if (xsProduct != 0.)
+  {
+    sigma = QuadInterpolator(  valueE11, valueE12, 
+    			       valueE21, valueE22, 
+			       xs11, xs12, 
+			       xs21, xs22, 
+			       valueT1, valueT2, 
+			       k, theta );
+  }
+
+  return sigma;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+
+G4double G4ecpssrLiCrossSection::FunctionFL2(G4double k, G4double theta) 							  
+{
+
+  G4double sigma = 0.;
+  G4double valueT1 = 0;
+  G4double valueT2 = 0;
+  G4double valueE21 = 0;
+  G4double valueE22 = 0;
+  G4double valueE12 = 0;
+  G4double valueE11 = 0;
+  G4double xs11 = 0;   
+  G4double xs12 = 0; 
+  G4double xs21 = 0; 
+  G4double xs22 = 0; 
+
+  // PROTECTION TO ALLOW INTERPOLATION AT MINIMUM AND MAXIMUM EtaK/Theta2 values 
+  // (in particular for FK computation at 95 for high velocity formula)
+  
+  if (
+  theta==9.5e-2 ||
+  theta==9.5e-1 ||
+  theta==9.5e+00 ||
+  theta==9.5e+01
+  ) theta=theta-1e-12;
+
+  if (
+  theta==1.e-2 ||
+  theta==1.e-1 ||
+  theta==1.e+00 ||
+  theta==1.e+01
+  ) theta=theta+1e-12;
+
+  // END PROTECTION
+  
+  {
+    std::vector<double>::iterator t2 = std::upper_bound(dummyVec.begin(),dummyVec.end(), k);
+    std::vector<double>::iterator t1 = t2-1;
+ 
+    std::vector<double>::iterator e12 = std::upper_bound(aVecMap[(*t1)].begin(),aVecMap[(*t1)].end(), theta);
+    std::vector<double>::iterator e11 = e12-1; 
+	  
+    std::vector<double>::iterator e22 = std::upper_bound(aVecMap[(*t2)].begin(),aVecMap[(*t2)].end(), theta);
+    std::vector<double>::iterator e21 = e22-1;
+	  	
+    valueT1  =*t1;
+    valueT2  =*t2;
+    valueE21 =*e21;
+    valueE22 =*e22;
+    valueE12 =*e12;
+    valueE11 =*e11;
+
+    xs11 = FL2Data[valueT1][valueE11];
+    xs12 = FL2Data[valueT1][valueE12];
+    xs21 = FL2Data[valueT2][valueE21];
+    xs22 = FL2Data[valueT2][valueE22];
+
+
+}
+     
+  G4double xsProduct = xs11 * xs12 * xs21 * xs22;
+  
+  if (xs11==0 || xs12==0 ||xs21==0 ||xs22==0) return (0.);
+     
+  if (xsProduct != 0.)
+  {
+    sigma = QuadInterpolator(  valueE11, valueE12, 
+    			       valueE21, valueE22, 
+			       xs11, xs12, 
+			       xs21, xs22, 
+			       valueT1, valueT2, 
+			       k, theta );
+  }
+
+  return sigma;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ecpssrLiCrossSection::LinLogInterpolate(G4double e1, 
+						        G4double e2, 
+						        G4double e, 
+						        G4double xs1, 
+						        G4double xs2)
+{
+  G4double d1 = std::log(xs1);
+  G4double d2 = std::log(xs2);
+  G4double value = std::exp(d1 + (d2 - d1)*(e - e1)/ (e2 - e1));
+  return value;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ecpssrLiCrossSection::LogLogInterpolate(G4double e1, 
+						        G4double e2, 
+						        G4double e, 
+						        G4double xs1, 
+						        G4double xs2)
+{
+  G4double a = (std::log10(xs2)-std::log10(xs1)) / (std::log10(e2)-std::log10(e1));
+  G4double b = std::log10(xs2) - a*std::log10(e2);
+  G4double sigma = a*std::log10(e) + b;
+  G4double value = (std::pow(10.,sigma));
+  return value;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4ecpssrLiCrossSection::QuadInterpolator(G4double e11, G4double e12, 
+						       G4double e21, G4double e22, 
+						       G4double xs11, G4double xs12, 
+						       G4double xs21, G4double xs22, 
+						       G4double t1, G4double t2, 
+						       G4double t, G4double e)
+{
+// Log-Log
+/*
+  G4double interpolatedvalue1 = LogLogInterpolate(e11, e12, e, xs11, xs12);
+  G4double interpolatedvalue2 = LogLogInterpolate(e21, e22, e, xs21, xs22);
+  G4double value = LogLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
+*/
+
+// Lin-Log
+  G4double interpolatedvalue1 = LinLogInterpolate(e11, e12, e, xs11, xs12);
+  G4double interpolatedvalue2 = LinLogInterpolate(e21, e22, e, xs21, xs22);
+  G4double value = LinLogInterpolate(t1, t2, t, interpolatedvalue1, interpolatedvalue2);
+  return value;
+}
+
Index: trunk/source/processes/electromagnetic/lowenergy/src/G4empCrossSection.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4empCrossSection.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4empCrossSection.cc	(revision 1197)
@@ -0,0 +1,123 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//$Id: G4empCrossSection.cc,v 1.2 2009/06/25 15:52:08 mantero Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+//         
+//
+// History:
+// -----------
+//  29 Apr 2009   ALF   1st implementation
+//
+// -------------------------------------------------------------------
+// Class description:
+// empirical model for K and L Ionization CS for Protons and Alpha
+// Further documentation available from http://www.ge.infn.it/geant4/lowE
+// -------------------------------------------------------------------
+
+
+#include "globals.hh"
+#include "G4empCrossSection.hh"
+#include "G4Proton.hh"
+//#include "G4Alpha.hh"
+//#include <math.h>
+
+G4empCrossSection::G4empCrossSection()
+  :totalCS(0)
+{ 
+
+  paulShellK = new G4PaulKCrossSection();
+  orlicShellLi = new G4OrlicLiCrossSection();
+
+}
+
+G4empCrossSection::~G4empCrossSection()
+{ 
+
+  delete paulShellK;
+  delete orlicShellLi;
+
+}
+
+std::vector<G4double> G4empCrossSection::GetCrossSection(G4int Z,
+							     G4double incidentEnergy,
+							     G4double mass,
+							     G4double deltaEnergy,
+							     G4bool testFlag) const
+{
+
+  deltaEnergy = 0;
+  testFlag = 0;
+
+  std::vector<G4double> crossSections;
+
+  crossSections.push_back( paulShellK->CalculateKCrossSection(Z, mass, incidentEnergy) );
+  
+  G4Proton* aProtone = G4Proton::Proton();
+  
+  if (mass == aProtone->GetPDGMass() ) {
+    
+    crossSections.push_back( orlicShellLi->CalculateL1CrossSection(Z, incidentEnergy) );
+    crossSections.push_back( orlicShellLi->CalculateL2CrossSection(Z, incidentEnergy) );
+    crossSections.push_back( orlicShellLi->CalculateL3CrossSection(Z, incidentEnergy) );
+  }
+  
+  return crossSections;
+
+}
+
+
+
+
+std::vector<G4double> G4empCrossSection::Probabilities(G4int Z,
+							   G4double incidentEnergy,
+							   G4double mass,
+							   G4double deltaEnergy) const
+{
+  
+std::vector<G4double> crossSections = GetCrossSection(Z, incidentEnergy, mass, deltaEnergy);
+
+  for (size_t i=0; i<crossSections.size(); i++ ) {
+    
+    if (totalCS) {
+      crossSections[i] = crossSections[i]/totalCS;
+    }
+    
+  }
+
+  return crossSections;
+
+}
+
+
+void G4empCrossSection::SetTotalCS(G4double val){
+
+  totalCS = val;
+
+}
+
+
+
Index: trunk/source/processes/electromagnetic/lowenergy/src/G4teoCrossSection.cc
===================================================================
--- trunk/source/processes/electromagnetic/lowenergy/src/G4teoCrossSection.cc	(revision 1197)
+++ trunk/source/processes/electromagnetic/lowenergy/src/G4teoCrossSection.cc	(revision 1197)
@@ -0,0 +1,145 @@
+//
+// ********************************************************************
+// * License and Disclaimer                                           *
+// *                                                                  *
+// * The  Geant4 software  is  copyright of the Copyright Holders  of *
+// * the Geant4 Collaboration.  It is provided  under  the terms  and *
+// * conditions of the Geant4 Software License,  included in the file *
+// * LICENSE and available at  http://cern.ch/geant4/license .  These *
+// * include a list of copyright holders.                             *
+// *                                                                  *
+// * Neither the authors of this software system, nor their employing *
+// * institutes,nor the agencies providing financial support for this *
+// * work  make  any representation or  warranty, express or implied, *
+// * regarding  this  software system or assume any liability for its *
+// * use.  Please see the license in the file  LICENSE  and URL above *
+// * for the full disclaimer and the limitation of liability.         *
+// *                                                                  *
+// * This  code  implementation is the result of  the  scientific and *
+// * technical work of the GEANT4 collaboration.                      *
+// * By using,  copying,  modifying or  distributing the software (or *
+// * any work based  on the software)  you  agree  to acknowledge its *
+// * use  in  resulting  scientific  publications,  and indicate your *
+// * acceptance of all terms of the Geant4 Software license.          *
+// ********************************************************************
+//
+//$Id: G4teoCrossSection.cc,v 1.5 2009/11/11 09:14:53 mantero Exp $
+// GEANT4 tag $Name: geant4-09-03-cand-01 $
+//
+//         
+//
+// History:
+// -----------
+//  21 Apr 2009   ALF   1st implementation
+//  29 Apr 2009   ALF Updated Desing for Integration
+//  11 Nov 2009   ALF update and code cleaning for the Dec Release 
+//
+// -------------------------------------------------------------------
+// Class description:
+// Low Energy Electromagnetic Physics, Cross section, p ionisation, K shell
+// Further documentation available from http://www.ge.infn.it/geant4/lowE
+// -------------------------------------------------------------------
+
+
+#include "globals.hh"
+#include "G4teoCrossSection.hh"
+//#include "G4AtomicTransitionManager.hh"
+//#include "G4NistManager.hh"
+#include "G4Proton.hh"
+//#include "G4Alpha.hh"
+//#include <math.h>
+
+G4teoCrossSection::G4teoCrossSection(G4String shellModel = "ecpssr")
+  :totalCS(0)
+{ 
+
+  if (shellModel == "ecpssr") {
+
+    ecpssrShellK = new G4ecpssrKCrossSection();
+
+    ecpssrShellLi = new G4ecpssrLiCrossSection();
+  }
+  else {
+
+    G4cout << "G4teoCrossSection: sorry, only ECPSSR model available at the moment" << G4endl;
+
+    ecpssrShellK = new G4ecpssrKCrossSection();
+
+    ecpssrShellLi = new G4ecpssrLiCrossSection();
+
+  }
+
+
+}
+
+G4teoCrossSection::~G4teoCrossSection()
+{ 
+
+  delete ecpssrShellK;
+  delete ecpssrShellLi;
+
+}
+
+std::vector<G4double> G4teoCrossSection::GetCrossSection(G4int Z,
+							     G4double incidentEnergy,
+							     G4double mass,
+							     G4double deltaEnergy,
+							     G4bool testFlag) const
+{
+
+  deltaEnergy = 0;
+  testFlag = 0;
+
+
+  std::vector<G4double> crossSections;
+
+  crossSections.push_back( ecpssrShellK->CalculateCrossSection(Z, mass, incidentEnergy) );
+  
+  //  G4Proton* aProtone = G4Proton::Proton();
+  
+  //  if (mass == aProtone->GetPDGMass() ) {
+    
+
+  //  }
+  
+    crossSections.push_back( ecpssrShellLi->CalculateL1CrossSection(Z, mass, incidentEnergy) );
+    crossSections.push_back( ecpssrShellLi->CalculateL2CrossSection(Z, mass, incidentEnergy) );
+    crossSections.push_back( ecpssrShellLi->CalculateL3CrossSection(Z, mass, incidentEnergy) );
+
+
+  return crossSections;
+
+}
+
+
+
+
+std::vector<G4double> G4teoCrossSection::Probabilities(G4int Z,
+							   G4double incidentEnergy,
+							   G4double mass,
+							   G4double deltaEnergy) const
+{
+  
+std::vector<G4double> crossSections = GetCrossSection(Z, incidentEnergy, mass, deltaEnergy);
+
+  for (size_t i=0; i<crossSections.size(); i++ ) {
+    
+    if (totalCS) {
+      crossSections[i] = crossSections[i]/totalCS;
+    }
+    
+  }
+
+  return crossSections;
+
+}
+
+
+void G4teoCrossSection::SetTotalCS(G4double val){
+
+  totalCS = val;
+
+}
+
+
+
