Changeset 961 for trunk/source/processes/electromagnetic/utils/src

Timestamp:

Apr 6, 2009, 12:21:12 PM (15 years ago)

Author:

garnier

Message:

update processes

Location:

trunk/source/processes/electromagnetic/utils/src

Files:

• G4DummyModel.cc (modified) (1 diff)
• G4EmCalculator.cc (modified) (9 diffs)
• G4EmCorrections.cc (modified) (15 diffs)
• G4EmModelManager.cc (modified) (29 diffs)
• G4EmMultiModel.cc (modified) (1 diff)
• G4EmProcessOptions.cc (modified) (5 diffs)
• G4EnergyLossMessenger.cc (modified) (19 diffs)
• G4EnergyLossTables.cc (modified) (2 diffs)
• G4LossTableBuilder.cc (modified) (9 diffs)
• G4LossTableManager.cc (modified) (16 diffs)
• G4VEmFluctuationModel.cc (modified) (3 diffs)
• G4VEmModel.cc (modified) (6 diffs)
• G4VEmProcess.cc (modified) (12 diffs)
• G4VEnergyLoss.cc (modified) (1 diff)
• G4VEnergyLossProcess.cc (modified) (50 diffs)
• G4VMultipleScattering.cc (modified) (14 diffs)
• G4ionEffectiveCharge.cc (modified) (10 diffs)

trunk/source/processes/electromagnetic/utils/src/G4DummyModel.cc

@@ -25 +25 @@
 //
 // $Id: G4DummyModel.cc,v 1.3 2007/05/22 17:31:58 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------

trunk/source/processes/electromagnetic/utils/src/G4EmCalculator.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmCalculator.cc,v 1.37 2007/08/16 15:55:42 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4EmCalculator.cc,v 1.46 2009/02/24 09:56:03 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -55 +55 @@
 //            10 keV to 1 keV (V. Ivanchenko)
 // 15.03.2007 Add ComputeEnergyCutFromRangeCut methods (V.Ivanchenko)
+// 21.04.2008 Updated computations for ions (V.Ivanchenko)
 //
 // Class Description:
@@ -126 +127 @@
   if(couple && UpdateParticle(p, kinEnergy) ) {
     res = manager->GetDEDX(p, kinEnergy, couple);
+    if(isIon) {
+      G4double eth = 2.0*MeV/massRatio;
+      if(kinEnergy > eth) {
+        G4double x1 = corr->ComputeIonCorrections(p,mat,kinEnergy);
+        G4double x2 = corr->ComputeIonCorrections(p,mat,eth);
+        res += x1 - x2*eth/kinEnergy;
+        /*      
+        G4cout << "### GetDEDX: E= " << kinEnergy << " res= " << res 
+               << " x1= " << x1 << " x2= " << x2 
+               << " del= " << x1 - x2*eth/kinEnergy << G4endl;;
+        */
+      } 
+    }
+
     if(verbose>0) {
       G4cout << "G4EmCalculator::GetDEDX: E(MeV)= " << kinEnergy/MeV
@@ -406 +421 @@
         res = currentModel->ComputeDEDXPerVolume(
               mat, baseParticle, escaled, cut) * chargeSquare;
-        if(verbose > 1)
+        if(verbose > 1) {
           G4cout <<  baseParticle->GetParticleName()
                  << " Escaled(MeV)= " << escaled;
+        }
       } else {
         res = currentModel->ComputeDEDXPerVolume(mat, p, kinEnergy, cut);
         if(verbose > 1) G4cout <<  " no basePart E(MeV)= " << kinEnergy;
       }
-      if(verbose > 1)
+      if(verbose > 1) {
         G4cout << " DEDX(MeV/mm)= " << res*mm/MeV
                << " DEDX(MeV*cm^2/g)= "
                << res*gram/(MeV*cm2*mat->GetDensity())
                << G4endl;
-
-      if(isIon) {
-        if(currentModel->HighEnergyLimit() > 100.*MeV)
-          res += corr->HighOrderCorrections(p,mat,kinEnergy);
-        else
-          res *= corr->EffectiveChargeCorrection(p,mat,kinEnergy);
-        if(verbose > 1)
-          G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
-                 << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
-                 << G4endl;
-      }
-
-      
+      }
+
       // emulate boundary region for different parameterisations
       G4double eth = currentModel->LowEnergyLimit();
@@ -447 +452 @@
           res0 = loweModel->ComputeDEDXPerVolume(mat, p, eth, cut);
         }
-        if(verbose > 1)
+        if(verbose > 1) {
           G4cout << "At boundary energy(MeV)= " << eth/MeV
                  << " DEDX(MeV/mm)= " << res1*mm/MeV
                  << G4endl;
-        if(isIon) res1 += corr->HighOrderCorrections(p,mat,eth/massRatio);
-        //G4cout << "eth= " << eth << " escaled= " << escaled
-        //  << " res0= " << res0 << " res1= "
-        //       << res1 <<  "  q2= " << chargeSquare << G4endl;
+        }
+        /*
+        G4cout << "eth= " << eth << " escaled= " << escaled
+               << " res0= " << res0 << " res1= "
+               << res1 <<  "  q2= " << chargeSquare << G4endl;
+        */
         res *= (1.0 + (res0/res1 - 1.0)*eth/escaled);
+
+        if(isIon) {
+          G4double ethscaled = eth/massRatio;
+          if(kinEnergy > ethscaled) {
+            G4double x1 = corr->ComputeIonCorrections(p,mat,kinEnergy);
+            G4double x2 = corr->ComputeIonCorrections(p,mat,ethscaled);
+            res += x1 - x2*ethscaled/kinEnergy;
+          } 
+        
+          if(verbose > 1) {
+            G4cout << "After Corrections: DEDX(MeV/mm)= " << res*mm/MeV
+                   << " DEDX(MeV*cm^2/g)= " << res*gram/(MeV*cm2*mat->GetDensity())
+                   << G4endl;
+          }
+        }
       }
       
@@ -507 +529 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, const G4ParticleDefinition* part,
-                                          const G4Material* mat, G4double cut)
+G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, 
+                                          const G4ParticleDefinition* part,
+                                          const G4Material* mat, 
+                                          G4double cut)
 {
   G4double dedx = ComputeElectronicDEDX(kinEnergy,part,mat,cut);
@@ -517 +541 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, const G4String& part,
-                                          const G4String& mat, G4double cut)
+G4double G4EmCalculator::ComputeTotalDEDX(G4double kinEnergy, 
+                                          const G4String& part,
+                                          const G4String& mat, 
+                                          G4double cut)
 {
   return ComputeTotalDEDX(kinEnergy,FindParticle(part),FindMaterial(mat),cut);
@@ -729 +755 @@
     if(currentProcess) currentProcessName = currentProcess->GetProcessName();
 
-    if(p->GetParticleType() == "nucleus" &&
-       currentParticleName  != "deuteron" && currentParticleName != "triton") {
+    if(p->GetParticleType() == "nucleus" 
+       && currentParticleName != "deuteron"  
+       && currentParticleName != "triton"
+       && currentParticleName != "alpha+"
+       && currentParticleName != "helium"
+       && currentParticleName != "hydrogen"
+      ) {
       baseParticle = theGenericIon;
       massRatio = baseParticle->GetPDGMass()/p->GetPDGMass();
@@ -755 +786 @@
   if(isIon) {
     chargeSquare =
-     ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy);
-    if(currentProcess)
+     ionEffCharge->EffectiveChargeSquareRatio(p, currentMaterial, kinEnergy)
+      * corr->EffectiveChargeCorrection(p,currentMaterial,kinEnergy);
+    if(currentProcess) {
       currentProcess->SetDynamicMassCharge(massRatio,chargeSquare);
-    // G4cout << "massR= " << massRatio << "   q2= " << chargeSquare << G4endl;
+      //G4cout << "NewP: massR= " << massRatio << "   q2= " << chargeSquare << G4endl;
+    }
   }
   return true;

trunk/source/processes/electromagnetic/utils/src/G4EmCorrections.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmCorrections.cc,v 1.23.2.1 2008/04/22 15:28:13 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4EmCorrections.cc,v 1.51 2008/12/18 13:01:44 gunter Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -44 +44 @@
 // 08.05.2007 V.Ivanchenko Use G4IonTable for ion mass instead of NistTable to avoid
 //                         division by zero
+// 29.02.2008 V.Ivanchenko use expantions for log and power function
+// 21.04.2008 Updated computations for ions (V.Ivanchenko)
+// 20.05.2008 Removed Finite Size correction (V.Ivanchenko)
 //
 //
@@ -56 +59 @@
 #include "G4EmCorrections.hh"
 #include "Randomize.hh"
-#include "G4NistManager.hh"
 #include "G4ParticleTable.hh"
 #include "G4IonTable.hh"
 #include "G4VEmModel.hh"
 #include "G4Proton.hh"
+#include "G4GenericIon.hh"
 #include "G4LPhysicsFreeVector.hh"
+#include "G4PhysicsLogVector.hh"
+#include "G4ProductionCutsTable.hh"
+#include "G4MaterialCutsCouple.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -67 +73 @@
 G4EmCorrections::G4EmCorrections()
 {
-  Initialise();
   particle   = 0;
   curParticle= 0;
@@ -74 +79 @@
   curVector  = 0;
   kinEnergy  = 0.0;
-  ionModel   = 0;
+  ionLEModel = 0;
+  ionHEModel = 0;
   nIons      = 0;
   verbose    = 1;
+  ncouples   = 0;
   massFactor = 1.0;
+  eth        = 2.0*MeV;
+  nbinCorr   = 20;
+  eCorrMin   = 25.*keV;
+  eCorrMax   = 250.*MeV;
   nist = G4NistManager::Instance();
   ionTable = G4ParticleTable::GetParticleTable()->GetIonTable();
+  Initialise();
 }
 
@@ -93 +105 @@
 G4double G4EmCorrections::HighOrderCorrections(const G4ParticleDefinition* p,
                                                const G4Material* mat,
-                                               G4double e)
+                                               G4double e, G4double)
 {
 // . Z^3 Barkas effect in the stopping power of matter for charged particles
@@ -108 +120 @@
   G4double Bloch  = BlochCorrection (p, mat, e);
   G4double Mott   = MottCorrection (p, mat, e);
-  G4double FSize  = FiniteSizeCorrection (p, mat, e);
-
-  G4double sum = (2.0*(Barkas + Bloch) + FSize + Mott);
+
+  G4double sum = (2.0*(Barkas + Bloch) + Mott);
 
   if(verbose > 1)
     G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
-           << " Bloch= " << Bloch << " Mott= " << Mott << " Fsize= " << FSize
+           << " Bloch= " << Bloch << " Mott= " << Mott 
            << " Sum= " << sum << G4endl; 
 
   sum *= material->GetElectronDensity() * q2 *  twopi_mc2_rcl2 /beta2;
+  return sum;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4EmCorrections::IonBarkasCorrection(const G4ParticleDefinition* p,
+                                              const G4Material* mat,
+                                              G4double e)
+{
+// . Z^3 Barkas effect in the stopping power of matter for charged particles
+//   J.C Ashley and R.H.Ritchie
+//   Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
+//   and ICRU49 report
+//   valid for kineticEnergy < 0.5 MeV
+
+  SetupKinematics(p, mat, e);
+  G4double res = 0.0;
+  if(tau > 0.0) 
+    res = 2.0*BarkasCorrection(p, mat, e)*
+      material->GetElectronDensity() * q2 *  twopi_mc2_rcl2 /beta2;
+  return res;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4EmCorrections::ComputeIonCorrections(const G4ParticleDefinition* p,
+                                                const G4Material* mat,
+                                                G4double e)
+{
+// . Z^3 Barkas effect in the stopping power of matter for charged particles
+//   J.C Ashley and R.H.Ritchie
+//   Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
+//   and ICRU49 report
+//   valid for kineticEnergy < 0.5 MeV
+//   Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
+  SetupKinematics(p, mat, e);
+  if(tau <= 0.0) return 0.0;
+
+  G4double Barkas = BarkasCorrection (p, mat, e);
+  G4double Bloch  = BlochCorrection (p, mat, e);
+  G4double Mott   = MottCorrection (p, mat, e);
+
+  G4double sum = 2.0*(Barkas*(charge - 1.0)/charge + Bloch) + Mott;
+
+  if(verbose > 1) {
+    G4cout << "EmCorrections: E(MeV)= " << e/MeV << " Barkas= " << Barkas
+           << " Bloch= " << Bloch << " Mott= " << Mott 
+           << " Sum= " << sum << G4endl; 
+  }
+  sum *= material->GetElectronDensity() * q2 *  twopi_mc2_rcl2 /beta2;
+
+  if(verbose > 1) G4cout << " Sum= " << sum << G4endl; 
+  return sum;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4EmCorrections::IonHighOrderCorrections(const G4ParticleDefinition* p,
+                                                  const G4MaterialCutsCouple* couple,
+                                                  G4double e)
+{
+// . Z^3 Barkas effect in the stopping power of matter for charged particles
+//   J.C Ashley and R.H.Ritchie
+//   Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
+//   and ICRU49 report
+//   valid for kineticEnergy < 0.5 MeV
+//   Other corrections from S.P.Ahlen Rev. Mod. Phys., Vol 52, No1, 1980
+
+  G4double sum = 0.0;
+
+  if(ionHEModel) {
+    G4int Z = G4int(p->GetPDGCharge()/eplus + 0.5);
+    if(Z >= 100)   Z = 99;
+    else if(Z < 1) Z = 1;
+
+    // fill vector
+    if(thcorr[Z].size() == 0) {
+      thcorr[Z].resize(ncouples);
+      G4double ethscaled = eth*p->GetPDGMass()/proton_mass_c2;
+
+      for(size_t i=0; i<ncouples; i++) {
+        (thcorr[Z])[i] = ethscaled*ComputeIonCorrections(p, currmat[i], ethscaled);
+        //G4cout << i << ". ethscaled= " << ethscaled 
+        //<< " corr= " << (thcorr[Z])[i]/ethscaled << G4endl;
+      } 
+    } 
+    G4double rest = (thcorr[Z])[couple->GetIndex()];
+
+    sum = ComputeIonCorrections(p,couple->GetMaterial(),e) - rest/e;
+
+    if(verbose > 1) G4cout << " Sum= " << sum << " dSum= " << rest/e << G4endl; 
+  }
   return sum;
 }
@@ -226 +329 @@
     G4int ieta = Index(y, Eta, nEtaK);
     corr = Value2(x, y, TheK[itet], TheK[itet+1], Eta[ieta], Eta[ieta+1],
-                  CK[itet][ieta], CK[itet+1][ieta], CK[itet][ieta+1], CK[itet+1][ieta+1]);
+                  CK[itet][ieta], CK[itet+1][ieta], 
+                  CK[itet][ieta+1], CK[itet+1][ieta+1]);
     //G4cout << "   x= " <<x<<" y= "<<y<<" tet= " <<TheK[itet]
     //<<" "<< TheK[itet+1]<<" eta= "<< Eta[ieta]<<" "<< Eta[ieta+1]
@@ -439 +543 @@
     else {
       G4int iw = Index(W, engBarkas, 47);
-      val = Value(W, engBarkas[iw], engBarkas[iw+1], corBarkas[iw], corBarkas[iw+1]);
+      val = Value(W, engBarkas[iw], engBarkas[iw+1], 
+                  corBarkas[iw], corBarkas[iw+1]);
     }
     //    G4cout << "i= " << i << " b= " << b << " W= " << W 
@@ -479 +584 @@
 {
   SetupKinematics(p, mat, e);
-  G4double mterm = pi*fine_structure_const*beta*charge;
-  /*
-  G4double mterm = 0.0; 
-  if(beta > 0.0) {
-
-    // Estimation of mean square root of the ionisation potential 
-    G4double Zeff = 0.0;
-    G4double norm = 0.0;
-
-    for (G4int i = 0; i<numberOfElements; i++) {
-      G4double Z = (*theElementVector)[i]->GetZ();
-      Zeff += Z*atomDensity[i];
-      norm += atomDensity[i];
-    }
-    Zeff *= (2.0/norm);
-    G4double ze1 = std::log(Zeff);
-    G4double eexc= material->GetIonisation()->GetMeanExcitationEnergy()*ze1*ze1/Zeff;
-
-    G4double invbeta = 1.0/beta;
-    G4double invbeta2= invbeta*invbeta;
-    G4double za  = charge*fine_structure_const;
-    G4double za2 = za*za;
-    G4double za3 = za2*za;
-    G4double x   = za*invbeta;
-    G4double cosx;
-    if(x < COSEB[13]) {
-      G4int i = Index(x,COSEB,14);
-      cosx    = Value(x,COSEB[i], COSEB[i+1],COSXI[i],COSXI[i+1]);
-    } else {
-      cosx = COSXI[13]*COSEB[13]/x;
-    }
-
-    mterm =
-      za*beta*(1.725 + pi*cosx*(0.52 - 2.0*std::sqrt(eexc/(2.0*electron_mass_c2*bg2))));
-      + za2*(3.246 - 0.451*beta2)
-      + za3*(1.522*beta + 0.987*invbeta)
-      + za2*za2*(4.569 - 0.494*beta2 - 2.696*invbeta2)
-      + za3*za2*(1.254*beta + 0.222*invbeta - 1.17*invbeta*invbeta2);
-  }
-*/  
+  G4double mterm = CLHEP::pi*fine_structure_const*beta*charge;
   return mterm;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4EmCorrections::FiniteSizeCorrection(const G4ParticleDefinition* p,
-                                               const G4Material* mat,
-                                               G4double e)
-  // Finite size corrections are parameterized according to
-  // J.D.Jackson Phys. Rev. D59 (1998) 017301 
-{
-  if(p) return 0.0;
-  SetupKinematics(p, mat, e);
-  G4double term = 0.0;
-  G4double numlim = 0.2;
-
-  //Leptons
-  if(p->GetLeptonNumber() != 0) {
-    G4double x =  tmax/(e + mass);
-    term = x*x;
-
-    // Pions and Kaons
-  } else if(p->GetPDGSpin() == 0.0 && q2 < 1.5) {
-    G4double xpi = 0.736*GeV;
-    G4double x   = 2.0*electron_mass_c2*tmax0/(xpi*xpi);
-    if(x <= numlim) term = -x*(1.0 - 0.5*x);
-    else            term = -std::log(1.0 + x);
-
-    // Protons and baryons
-  } else if(q2 < 1.5) {
-    G4double xp = 0.8426*GeV;
-    G4double x  = 2.0*electron_mass_c2*tmax0/(xp*xp);
-    G4double ksi2 = 2.79285*2.79285;
-    if(x <= numlim) term = -x*((1.0 + 5.0*x/6.0)/((1.0 + x)*(1 + x)) + 0.5*x);
-    else term = -x*(1.0 + 5.0*x/6.0)/((1.0 + x)*(1 + x)) - std::log(1.0 + x);
-    G4double b  = xp*0.5/mass;
-    G4double c  = xp*mass/(electron_mass_c2*(mass + e));
-    G4double lb = b*b;
-    G4double lb2= lb*lb;
-    G4double nu = 0.5*c*c;
-    G4double x1 = 1.0 + x;
-    G4double x2 = x1*x1;
-    G4double l1 = 1.0 - lb;
-    G4double l2 = l1*l1;
-    G4double lx = 1.0 + lb*x;
-    G4double ia = lb2*(lx*std::log(lx/x1)/x + l1 -
-                       0.5*x*l2/(lb*x1) + 
-                       x*(3.0 + 2.0*x)*l2*l1/(6.0*x2*lb2))/(l2*l2);
-    G4double ib = x*x*(3.0 + x)/(6.0*x2*x1); 
-    term += lb*((ksi2 - 1.0)*ia + nu*ksi2*ib);
-    // G4cout << "Proton F= " << term << " ia= " << ia << " ib= " << ib << " lb= " << lb<< G4endl;
-    
-    //ions
-  } else {
-    G4double xp = 0.8426*GeV/std::pow(mass/proton_mass_c2,-0.33333333);
-    G4double x  = 2.0*electron_mass_c2*tmax0/(xp*xp);
-    if(x <= numlim) term = -x*(1.0 - 0.5*x);
-    else            term = -std::log(1.0 + x);
-    //G4cout << "Ion F= " << term << G4endl;
-  }
-
-  return term;
 }
 
@@ -629 +634 @@
   if (er >= ed[0])       nloss = a[0];
   else {
+    // the table is inverse in energy
     for (G4int i=102; i>=0; i--)
     {
@@ -657 +663 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4ionEffectiveCharge* G4EmCorrections::GetIonEffectiveCharge(G4VEmModel* m)
-{
-  if(m) ionModel = m;
-  return &effCharge;
-}
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4int G4EmCorrections::GetNumberOfStoppingVectors()
-{
-  return nIons;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 G4double G4EmCorrections::EffectiveChargeCorrection(const G4ParticleDefinition* p,
                                                     const G4Material* mat,
@@ -676 +668 @@
 {
   G4double factor = 1.0;
-  if(p->GetPDGCharge() <= 2.5*eplus) return factor;
-  if(verbose > 1)
-    G4cout << "EffectiveChargeCorrection: " << p->GetParticleName() << " in " << mat->GetName()
+  if(p->GetPDGCharge() <= 2.5*eplus || nIons <= 0) return factor;
+  /*
+  if(verbose > 1) {
+    G4cout << "EffectiveChargeCorrection: " << p->GetParticleName() 
+           << " in " << mat->GetName()
            << " ekin(MeV)= " << ekin/MeV << G4endl;
+  }
+  */
   if(p != curParticle || mat != curMaterial) {
     curParticle = p;
-    curMaterial = 0;
+    curMaterial = mat;
     curVector   = 0;
-    G4int Z = p->GetAtomicNumber();
-    G4int A = p->GetAtomicMass();
-    if(verbose > 1) G4cout << "Zion= " << Z << " Aion= " << A << G4endl;
-    massFactor = proton_mass_c2/ionTable->GetIonMass(Z,A);
-    idx = 0;
-    for(; idx<nIons; idx++) {
-      if(Z == Zion[idx] && A == Aion[idx]) {
-        if(materialList[idx] == mat) {
-          curMaterial = mat;
-          curVector = stopData[idx];
-          break;
-        } else if(materialList[idx] == 0) {
-          if(materialName[idx] == mat->GetName())
-            curVector = InitialiseMaterial(mat);
-        }
+    currentZ = p->GetAtomicNumber();
+    if(verbose > 1) {
+      G4cout << "G4EmCorrections::EffectiveChargeCorrection: Zion= " 
+             << currentZ << " Aion= " << p->GetPDGMass()/amu_c2 << G4endl;
+    }
+    massFactor = proton_mass_c2/p->GetPDGMass();
+    idx = -1;
+    G4int dz = 1000;
+
+    for(G4int i=0; i<nIons; i++) {
+      if(materialList[i] == mat) {
+        G4int delz = currentZ - Zion[i];
+        if(delz < 0) delz = -delz;
+        if(delz < dz) {
+          idx = i;
+          dz = delz;
+          if(0 == delz) break;
+        }
       }
+    }
+    //    G4cout << " idx= " << idx << " dz= " << dz << G4endl;
+    if(idx > 0) {
+      if(!ionList[idx]) BuildCorrectionVector(); 
+      if(ionList[idx])  curVector = stopData[idx];
     }
   }
@@ -705 +709 @@
     G4bool b;
     factor = curVector->GetValue(ekin*massFactor,b);
-  }
-  if(verbose > 1) G4cout << "  factor= " << factor << G4endl;
+    if(verbose > 1) {
+      G4cout << "E= " << ekin << " factor= " << factor << " massfactor= " 
+             << massFactor << G4endl;
+    }
+  }
   return factor;
 }
@@ -714 +721 @@
 void G4EmCorrections::AddStoppingData(G4int Z, G4int A,
                                       const G4String& mname,
-                                      G4PhysicsVector& dVector)
-{
-  idx = 0;
-  for(; idx<nIons; idx++) {
-    if(Z == Zion[idx] && A == Aion[idx] && mname == materialName[idx])
-      break;
-  }
-  if(idx == nIons) {
+                                      G4PhysicsVector* dVector)
+{
+  G4int i = 0;
+  for(; i<nIons; i++) {
+    if(Z == Zion[i] && A == Aion[i] && mname == materialName[i]) break;
+  }
+  if(i == nIons) {
     Zion.push_back(Z);
     Aion.push_back(A);
     materialName.push_back(mname);
     materialList.push_back(0);
-    stopData.push_back(0);
+    ionList.push_back(0);
+    stopData.push_back(dVector);
     nIons++;
-  } else {
-    if(stopData[idx]) delete stopData[idx];
-  }
-  size_t nbins = dVector.GetVectorLength();
-  size_t n = 0;
-  for(; n<nbins; n++) {
-    if(dVector.GetLowEdgeEnergy(n) > 2.0*MeV) break;
-  }
-  if(n < nbins) nbins = n + 1;
-  G4LPhysicsFreeVector* v =
-    new G4LPhysicsFreeVector(nbins,
-                             dVector.GetLowEdgeEnergy(0),
-                             dVector.GetLowEdgeEnergy(nbins-1));
+    if(verbose>1) {
+      G4cout << "AddStoppingData Z= " << Z << " A= " << A << " " << mname
+             << "  idx= " << i << G4endl;
+    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4EmCorrections::BuildCorrectionVector()
+{
+  if(!ionLEModel || !ionHEModel) {
+    return;
+  }
+
+  const G4ParticleDefinition* ion = curParticle;
+  G4int Z = Zion[idx];
+  if(currentZ != Z) {
+    ion = G4ParticleTable::GetParticleTable()->FindIon(Z, Aion[idx], 0, Z);
+  }
+  //G4cout << "BuildCorrectionVector: idx= " << idx << " Z= " << Z 
+  //     << " curZ= " << currentZ << G4endl;
+
+  // G4double A = nist->GetAtomicMassAmu(Z);
+  G4double A = G4double(ion->GetBaryonNumber());
+  G4PhysicsVector* v = stopData[idx];
+    
+  const G4ParticleDefinition* p = G4GenericIon::GenericIon();
+  G4double massRatio = proton_mass_c2/ion->GetPDGMass();
+
+  if(verbose>1) {
+    G4cout << "BuildCorrectionVector: Stopping for "
+           << curParticle->GetParticleName() << " in " 
+           << materialName[idx] << " Ion Z= " << Z << " A= " << A
+           << " massRatio= " << massRatio << G4endl;
+  }
   G4bool b;
-  for(size_t i=0; i<nbins; i++) {
-    G4double e = dVector.GetLowEdgeEnergy(i);
-    G4double dedx = dVector.GetValue(e, b);
-    v->PutValues(i, e, dedx);
-  }
-  //  G4cout << *v << G4endl;
-  stopData[idx] = v;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4PhysicsVector* G4EmCorrections::InitialiseMaterial(const G4Material* mat)
-{
-  G4PhysicsVector* v = 0;
-  const G4Material* m = nist->FindOrBuildMaterial(materialName[idx],false);
-  if(m) {
-    materialList[idx] = m;
-    curMaterial = mat;
-    v = stopData[idx];
-    size_t nbins = v->GetVectorLength();
-    const G4ParticleDefinition* p = G4Proton::Proton();
-    if(verbose>1) G4cout << "G4ionIonisation::InitialiseMaterial Stooping data for "
-                         << materialName[idx] << G4endl;
-    G4bool b;
-    for(size_t i=0; i<nbins; i++) {
-      G4double e = v->GetLowEdgeEnergy(i);
-      G4double dedx = v->GetValue(e, b);
-      G4double dedx1= ionModel->ComputeDEDXPerVolume(mat, p, e, e)*
-        effCharge.EffectiveChargeSquareRatio(curParticle,mat,e/massFactor);
-      v->PutValue(i, dedx/dedx1);
-      if(verbose>1) G4cout << "  E(meV)= " << e/MeV << "   Correction= " << dedx/dedx1
-                           << "   "  << dedx << " " << dedx1 << G4endl;
-    }
-  }
-  return v;
+
+  G4PhysicsLogVector* vv = 
+    new G4PhysicsLogVector(eCorrMin,eCorrMax,nbinCorr);
+  vv->SetSpline(true);
+  G4double e, eion, dedx, dedx1;
+  G4double eth0 = v->GetLowEdgeEnergy(0);
+  G4double escal = eth/massRatio;
+  G4double qe = 
+    effCharge.EffectiveChargeSquareRatio(ion, curMaterial, escal); 
+  G4double dedxt = 
+    ionLEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe;
+  G4double dedx1t = 
+    ionHEModel->ComputeDEDXPerVolume(curMaterial, p, eth, eth)*qe 
+    + ComputeIonCorrections(curParticle, curMaterial, escal);
+  G4double rest = escal*(dedxt - dedx1t);
+  //G4cout << "Escal(MeV)= "<<escal<<" dedxt0= " <<dedxt 
+  //     << " dedxt1= " << dedx1t << G4endl;   
+
+  for(G4int i=0; i<nbinCorr; i++) {
+    e = vv->GetLowEdgeEnergy(i);
+    escal = e/massRatio;
+    eion  = escal/A;
+    if(eion <= eth0) {
+      dedx = v->GetValue(eth0, b)*std::sqrt(eion/eth0);
+    } else {
+      dedx = v->GetValue(eion, b);
+    }
+    qe = effCharge.EffectiveChargeSquareRatio(curParticle,curMaterial,escal); 
+    if(e <= eth) {
+      dedx1 = ionLEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe;
+    } else {
+      dedx1 = ionHEModel->ComputeDEDXPerVolume(curMaterial, p, e, e)*qe +
+        ComputeIonCorrections(curParticle, curMaterial, escal) + rest/escal;
+    }
+    vv->PutValue(i, dedx/dedx1);
+    if(verbose>1) {
+      G4cout << "  E(meV)= " << e/MeV << "   Correction= " << dedx/dedx1
+             << "   "  << dedx << " " << dedx1 
+             << "  massF= " << massFactor << G4endl;
+    }
+  }
+  delete v;
+  ionList[idx]  = ion;
+  stopData[idx] = vv;
+  if(verbose>1) G4cout << "End data set " << G4endl;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4EmCorrections::InitialiseForNewRun()
+{
+  G4ProductionCutsTable* tb = G4ProductionCutsTable::GetProductionCutsTable();
+  ncouples = tb->GetTableSize();
+  if(currmat.size() != ncouples) {
+    currmat.resize(ncouples);
+    size_t i;
+    for(i=0; i<100; i++) {thcorr[i].clear();}
+    for(i=0; i<ncouples; i++) {
+      currmat[i] = tb->GetMaterialCutsCouple(i)->GetMaterial();
+      G4String nam = currmat[i]->GetName();
+      for(G4int j=0; j<nIons; j++) {
+        if(nam == materialName[j]) { materialList[j] = currmat[i]; }
+      }
+    }
+  }
 }
 

trunk/source/processes/electromagnetic/utils/src/G4EmModelManager.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmModelManager.cc,v 1.40 2007/11/09 11:35:54 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4EmModelManager.cc,v 1.46 2008/10/13 14:56:56 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -78 +78 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& list, G4DataVector& lowE)
+G4RegionModels::G4RegionModels(G4int nMod, std::vector<G4int>& indx, 
+                               G4DataVector& lowE, const G4Region* reg)
 {
   nModelsForRegion      = nMod;
   theListOfModelIndexes = new G4int [nModelsForRegion];
-  lowKineticEnergy      = new G4double [nModelsForRegion];
+  lowKineticEnergy      = new G4double [nModelsForRegion+1];
   for (G4int i=0; i<nModelsForRegion; i++) {
-    theListOfModelIndexes[i] = list[i];
+    theListOfModelIndexes[i] = indx[i];
     lowKineticEnergy[i] = lowE[i];
   }
+  lowKineticEnergy[nModelsForRegion] = lowE[nModelsForRegion];
+  theRegion = reg;
 }
 
@@ -99 +102 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-#include "G4LossTableManager.hh"
 #include "G4Step.hh"
 #include "G4ParticleDefinition.hh"
@@ -107 +109 @@
 #include "G4MaterialCutsCouple.hh"
 #include "G4ProductionCutsTable.hh"
-#include "G4Region.hh"
 #include "G4RegionStore.hh"
 #include "G4Gamma.hh"
 #include "G4Positron.hh"
+#include "G4UnitsTable.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -128 +130 @@
   regions.clear();
   orderOfModels.clear();
-  upperEkin.clear();
   maxCutInRange    = 12.*cm;
   maxSubCutInRange = 0.7*mm;
@@ -139 +140 @@
 G4EmModelManager::~G4EmModelManager()
 {
-  G4int i,j;
   Clear();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4EmModelManager::Clear()
+{
   if(1 < verboseLevel) {
-    G4cout << "G4EmModelManager:: delete models ";
-    if(particle) G4cout << " for " << particle->GetParticleName();
-    G4cout << " nModels=" << nEmModels <<G4endl;
-  }
-
-  for(i = 0; i<nEmModels; i++) {
-    orderOfModels[i] = 1;
-  }
-  for(i = 0; i<nEmModels; i++) {
-    if (orderOfModels[i]) {
-      orderOfModels[i] = 0;
-      for(j = i+1; j<nEmModels; j++) {
-        if(models[i] == models[j]) orderOfModels[j] = 0;
-      }
-      G4String nam = models[i]->GetName();
-      if(nam != "PAI" && nam != "PAIModel" ) delete models[i];
-    }
-  }
-  for(i = 0; i<nEmModels; i++) {
-    orderOfModels[i] = 1;
-  }
-  for(i = 0; i<nEmModels; i++) {
-    if (orderOfModels[i]) {
-      orderOfModels[i] = 0;
-      for(j = i+1; j<nEmModels; j++) {
-        if(flucModels[i] == flucModels[j]) orderOfModels[j] = 0;
-      }
-      delete flucModels[i];
-    }
-  }
-  if(1 < verboseLevel) 
-    G4cout << "G4EmModelManager:: models are deleted!" << G4endl;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4EmModelManager::Clear()
-{
-  if(1 < verboseLevel) {
-    G4cout << "G4EmModelManager::Clear()";
-    if(particle) G4cout << " for " << particle->GetParticleName();
-    G4cout << G4endl;
+    G4cout << "G4EmModelManager::Clear()" << G4endl;
   }
 
   theCuts.clear();
   theSubCuts.clear();
-  upperEkin.clear();
   if(idxOfRegionModels) delete [] idxOfRegionModels;
   if(setOfRegionModels && nRegions) {
@@ -215 +179 @@
   orderOfModels.push_back(num);
   p->DefineForRegion(r);
-  if (nEmModels>0) {
-    G4int idx = nEmModels;
-    do {idx--;} while (idx && num < orderOfModels[idx]);
-    if (num >= orderOfModels[idx] && num <= orderOfModels[idx+1]) idx++;
-    if (idx < nEmModels) {
-      models[nEmModels] = models[idx];
-      flucModels[nEmModels] = flucModels[idx];
-      regions[nEmModels] = regions[idx];
-      orderOfModels[nEmModels] = orderOfModels[idx];
-      models[idx] = p;
-      flucModels[idx] = fm;
-      regions[idx] = r;
-      orderOfModels[idx] = num;
-    }
-  }
   nEmModels++;
 }
@@ -254 +203 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4VEmModel* G4EmModelManager::GetModel(G4int i)
+G4VEmModel* G4EmModelManager::GetModel(G4int i, G4bool ver)
 {
   G4VEmModel* m = 0;
-  if(i >= 0 && i < nEmModels) m = models[i];
-  else if(verboseLevel > 0) 
+  if(i >= 0 && i < nEmModels) {m = models[i];}
+  else if(verboseLevel > 0 && ver) { 
     G4cout << "G4EmModelManager::GetModel WARNING: "
            << "index " << i << " is wrong Nmodels= "
-           << nEmModels << G4endl;
+           << nEmModels;
+    if(particle) G4cout << " for " << particle->GetParticleName(); 
+    G4cout<< G4endl;
+  }
   return m;
 }
@@ -269 +221 @@
 const G4DataVector* G4EmModelManager::Initialise(const G4ParticleDefinition* p,
                                                  const G4ParticleDefinition* sp,
-                                                       G4double theMinSubRange,
-                                                       G4int val)
+                                                 G4double theMinSubRange,
+                                                 G4int val)
 {
   verboseLevel = val;
@@ -293 +245 @@
   // Identify the list of regions with different set of models
   nRegions = 1;
-  std::vector<const G4Region*> set;
-  set.push_back(world);
+  std::vector<const G4Region*> setr;
+  setr.push_back(world);
   G4bool isWorld = false;
 
@@ -306 +258 @@
       if (nRegions>1) {
         for (G4int j=1; j<nRegions; j++) {
-          if ( r == set[j] ) newRegion = false;
+          if ( r == setr[j] ) newRegion = false;
         }
       }
       if (newRegion) {
-        set.push_back(r);
+        setr.push_back(r);
         nRegions++;
       }
@@ -322 +274 @@
   idxOfRegionModels[numOfCouples] = 0;
   setOfRegionModels = new G4RegionModels*[nRegions];
-  upperEkin.resize(nEmModels);
+
+  std::vector<G4int>    modelAtRegion(nEmModels);
+  std::vector<G4int>    modelOrd(nEmModels);
+  G4DataVector          eLow(nEmModels);
+  G4DataVector          eHigh(nEmModels);
+  G4int nmax = nEmModels;
 
   // Order models for regions
   for (G4int reg=0; reg<nRegions; reg++) {
-    const G4Region* region = set[reg];
-
+    const G4Region* region = setr[reg];
     G4int n = 0;
-
-    std::vector<G4int>    modelAtRegion;
-    G4DataVector            eLow;
-    G4DataVector            eHigh;
-    modelAtRegion.clear();
-    eLow.clear();
-    eHigh.clear();
 
     if(isWorld || 0 < reg) {
@@ -346 +295 @@
           G4double tmin = model->LowEnergyLimit();
           G4double tmax = model->HighEnergyLimit();
-          if (n>0) tmin = std::max(tmin, eHigh[n-1]);
+          G4int  ord    = orderOfModels[ii];
+          G4bool push   = true;
+          G4bool insert = false;
+          G4int  idx    = n;
 
           if(1 < verboseLevel) {
@@ -355 +307 @@
                  << " tmin(MeV)= " << tmin/MeV
                  << "; tmax(MeV)= " << tmax/MeV
-                 << "; order= " << orderOfModels[ii]
+                 << "; order= " << ord
                  << G4endl;
           }
-
-          if (tmin < tmax) {
-            modelAtRegion.push_back(ii);
-            eLow.push_back(tmin);
-            eHigh.push_back(tmax);
-            upperEkin[ii] = tmax;
-            n++;
+        
+          if (n == 0) n++;
+          else {
+            tmin = std::min(tmin, eHigh[n-1]);
+            if(tmin >= tmax) push = false;
+            else {
+
+              // high energy model
+              if(tmin == eHigh[n-1] && tmax > eHigh[n-1]) n++;
+              else if (tmax > eHigh[n-1]) {
+
+                // compare order of models
+                for(G4int k = n-1; k>=0; k--) {
+       
+                  if (ord >= modelOrd[k]) {
+                    tmin = std::max(tmin, eHigh[k]);           
+                    if(k < n-1) n = k + 2;
+                    break;
+                  } else if (tmin > eLow[k]) { 
+                    eHigh[k] = tmin;
+                    n = k + 2;
+                    break;
+                  } else if (tmin == eLow[k]) {
+                    n = k + 1;
+                    break;
+                  }
+                }
+                if(tmin < eLow[0]) n = 1;
+                idx = n - 1;
+
+                // low energy model
+              } else {
+
+                tmax = std::max(tmax, eLow[0]);
+                insert = true;
+                push = false;
+                idx = 0;
+                if(tmax <= eLow[0]) tmax = eLow[0];
+                else {
+
+                  for(G4int k=0; k<n; k++) {
+       
+                    if (ord >= modelOrd[k]) {
+                      if(k == 0) {
+                        if(tmin < eLow[0]) tmax = eLow[0];
+                        else insert = false;
+                        break; 
+                      } else {
+                        insert = false;
+                        break;
+                      }
+                    } else if(tmax < eHigh[k]) {
+                      idx = k;
+                      if(k > 0) tmin = eLow[k];
+                      eLow[k] = tmax;               
+                      break; 
+                    } else if(tmax == eHigh[k]) {           
+                      insert = false;
+                      push = true;
+                      idx = k;
+                      if(k > 0) tmin = eLow[k];
+                      else tmin = std::min(tmin,eLow[0]);
+                      break;
+                    } else {
+                      modelAtRegion[k] = ii;
+                      modelOrd[k] = ord;
+                      if(k == 0) eLow[idx] = std::min(tmin,eLow[0]);
+                    }
+                  }
+                }
+                if(insert && idx < n) n++;
+                else insert = false;
+              }
+            }
+          }
+          if(n > nmax) {
+            nmax = n;
+            modelAtRegion.resize(nmax);
+            modelOrd.resize(nmax);
+            eLow.resize(nmax);
+            eHigh.resize(nmax);
+          }
+          if(insert) {
+            for(G4int k=n-2; k>=idx; k--) {            
+              modelAtRegion[k+1] = modelAtRegion[k];
+              modelOrd[k+1] = modelOrd[k];
+              eLow[k+1]  = eLow[k];
+              eHigh[k+1] = eHigh[k];
+            }
+          }
+          if (push || insert) {
+            modelAtRegion[idx] = ii;
+            modelOrd[idx] = ord;
+            eLow[idx]  = tmin;
+            eHigh[idx] = tmax;
           }
         }
@@ -374 +414 @@
       eLow.push_back(0.0);
       eHigh.push_back(DBL_MAX);
-      upperEkin.push_back(DBL_MAX);
     }
     eLow[0] = 0.0;
+    if(n >= nmax) eLow.resize(nmax+1);
+    eLow[n] = eHigh[n-1];
 
     if(1 < verboseLevel) {
@@ -382 +423 @@
       if (region) G4cout << region->GetName();
       G4cout << ">  Elow(MeV)= ";
-      for(G4int ii=0; ii<n; ii++) {G4cout << eLow[ii]/MeV << " ";}
+      for(G4int ii=0; ii<=n; ii++) {G4cout << eLow[ii]/MeV << " ";}
       G4cout << G4endl;
     }
-    G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow);
+    G4RegionModels* rm = new G4RegionModels(n, modelAtRegion, eLow, region);
     setOfRegionModels[reg] = rm;
   }
@@ -399 +440 @@
  
     G4int reg = nRegions;
-    do {reg--;} while (reg>0 && pcuts != (set[reg]->GetProductionCuts()));
+    do {reg--;} while (reg>0 && pcuts != (setr[reg]->GetProductionCuts()));
     idxOfRegionModels[i] = reg;
 
@@ -474 +515 @@
 {
 
-  // vectors to provide continues dE/dx
-  G4DataVector factor;
-  G4DataVector dedxLow;
-  G4DataVector dedxHigh;
-
   G4double e;
 
@@ -491 +527 @@
     G4cout << "G4EmModelManager::FillDEDXVector() for "
            << couple->GetMaterial()->GetName()
-           << "  Ecut(MeV)= " << cut
-           << "  Esubcut(MeV)= " << subcut
+           << "  cut(MeV)= " << cut
+           << "  subcut(MeV)= " << subcut
            << "  Type " << tType
            << "  for " << particle->GetParticleName()
@@ -501 +537 @@
   const G4RegionModels* regModels = setOfRegionModels[reg];
   G4int nmod = regModels->NumberOfModels();
-  factor.resize(nmod);
-  dedxLow.resize(nmod);
-  dedxHigh.resize(nmod);
+
+  // vectors to provide continues dE/dx
+  G4DataVector factor(nmod);
+  G4DataVector eLow(nmod+1);
+  G4DataVector dedxLow(nmod);
+  G4DataVector dedxHigh(nmod);
 
   if(1 < verboseLevel) {
@@ -512 +551 @@
   }
 
-
   // calculate factors to provide continuity of energy loss
+
+
   factor[0] = 1.0;
   G4int j;
@@ -519 +559 @@
   G4int totBinsLoss = aVector->GetVectorLength();
 
-  dedxLow[0]  = 0.0;
-
-  e = upperEkin[regModels->ModelIndex(0)];
+  dedxLow[0] = 0.0;
+  eLow[0]    = 0.0;
+
+  e = regModels->LowEdgeEnergy(1);
+  eLow[1]    = e;
   G4VEmModel* model = models[regModels->ModelIndex(0)]; 
   dedxHigh[0] = 0.0;
   if(model && cut > subcut) {
     dedxHigh[0] = model->ComputeDEDX(couple,particle,e,cut);
-    if(subcut > 0.0) 
+    if(subcut > 0.0) {
       dedxHigh[0] -= model->ComputeDEDX(couple,particle,e,subcut);
+    }
   }
   if(nmod > 1) {
     for(j=1; j<nmod; j++) {
 
-      e = upperEkin[regModels->ModelIndex(j-1)];
+      e = regModels->LowEdgeEnergy(j);
+      eLow[j]   = e;
       G4int idx = regModels->ModelIndex(j); 
 
       dedxLow[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
-      if(subcut > 0.0) 
+      if(subcut > 0.0) {
         dedxLow[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
+      }
       if(subcut == cut) dedxLow[j] = 0.0;
 
-      e = upperEkin[idx];
+      e = regModels->LowEdgeEnergy(j+1);
+      eLow[j+1] = e;
       dedxHigh[j] = models[idx]->ComputeDEDX(couple,particle,e,cut);
-      if(subcut > 0.0) 
+      if(subcut > 0.0) {
         dedxHigh[j] -= models[idx]->ComputeDEDX(couple,particle,e,subcut);
+      }
       if(subcut == cut) dedxHigh[j] = 0.0;
     }
+    if(1 < verboseLevel) {
+      G4cout << " model #0"  
+             << "  dedx(" << eLow[0] << ")=  " << dedxLow[0]
+             << "  dedx(" << eLow[1] << ")=  " << dedxHigh[0]
+             << G4endl;
+    }
 
     for(j=1; j<nmod; j++) {
-      if(dedxLow[j] > 0.0) factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0);
-      else                 factor[j] = 0.0;
+      if(dedxLow[j] > 0.0) {
+        factor[j] = (dedxHigh[j-1]/dedxLow[j] - 1.0)*eLow[j];
+      } else  factor[j] = 0.0;
+      if(1 < verboseLevel) {
+        G4cout << " model #" << j 
+               << "  dedx(" << eLow[j] << ")=  " << dedxLow[j]
+               << "  dedx(" << eLow[j+1] << ")=  " << dedxHigh[j]
+               << "  factor= " << factor[j]/eLow[j]
+               << G4endl;
+      }
     }
 
@@ -565 +626 @@
       // Choose a model of energy losses
     G4int k = 0;
-    if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
+    if (nmod > 1 && e > eLow[1]) {
       do {
         k++;
-        fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
-      } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
+        fac *= (1.0 + factor[k]/e);
+      } while (k+1 < nmod && e > eLow[k+1]);
     }
 
@@ -602 +663 @@
                                         G4EmTableType tType)
 {
-  // vectors to provide continues cross section
-  G4DataVector factor;
-  G4DataVector sigmaLow;
-  G4DataVector sigmaHigh;
-
   G4double e;
 
@@ -630 +686 @@
   const G4RegionModels* regModels = setOfRegionModels[reg];
   G4int nmod = regModels->NumberOfModels();
-  factor.resize(nmod);
-  sigmaLow.resize(nmod);
-  sigmaHigh.resize(nmod);
+
+  // vectors to provide continues dE/dx
+  G4DataVector factor(nmod);
+  G4DataVector eLow(nmod+1);
+  G4DataVector sigmaLow(nmod);
+  G4DataVector sigmaHigh(nmod);
 
   if(2 < verboseLevel) {
@@ -644 +703 @@
   G4int totBinsLambda = aVector->GetVectorLength();
 
-  sigmaLow[0]  = 0.0;
-
-  e = upperEkin[regModels->ModelIndex(0)];
+  sigmaLow[0] = 0.0;
+  eLow[0]     = 0.0;
+
+  e = regModels->LowEdgeEnergy(1);
+  eLow[1]     = e;
   G4VEmModel* model = models[regModels->ModelIndex(0)];
   sigmaHigh[0] = 0.0;
@@ -668 +729 @@
     for(j=1; j<nmod; j++) {
 
-      e  = upperEkin[regModels->ModelIndex(j-1)];
-      sigmaLow[j] = 
-        models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
-      e  = upperEkin[regModels->ModelIndex(j)];
-      sigmaHigh[j] = 
-        models[regModels->ModelIndex(j)]->CrossSection(couple,particle,e,cut,tmax);
+      e  = regModels->LowEdgeEnergy(j);
+      eLow[j]   = e;
+      G4int idx = regModels->ModelIndex(j); 
+
+      sigmaLow[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
+      e  = regModels->LowEdgeEnergy(j+1);
+      eLow[j+1] = e;
+      sigmaHigh[j] = models[idx]->CrossSection(couple,particle,e,cut,tmax);
+    }
+    if(1 < verboseLevel) {
+      G4cout << " model #0"  
+             << "  sigma(" << eLow[0] << ")=  " << sigmaLow[0]
+             << "  sigma(" << eLow[1] << ")=  " << sigmaHigh[0]
+             << G4endl;
+    }
+    for(j=1; j<nmod; j++) {
+      if(sigmaLow[j] > 0.0) {
+        factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0)*eLow[j];
+      } else  factor[j] = 0.0;
       if(1 < verboseLevel) {
-        G4cout << " model #" << j << "   eUp= " << e 
-               << "  sigmaUp= " << sigmaHigh[j]
-               << "  sigmaDown= " << sigmaLow[j]
+        G4cout << " model #" << j 
+               << "  sigma(" << eLow[j] << ")=  " << sigmaLow[j]
+               << "  sigma(" << eLow[j+1] << ")=  " << sigmaHigh[j]
+               << "  factor= " << factor[j]/eLow[j]
                << G4endl;
       }
-    }
-    for(j=1; j<nmod; j++) {
-      if(sigmaLow[j] > 0.0) factor[j] = (sigmaHigh[j-1]/sigmaLow[j] - 1.0);
-      else                  factor[j] = 0.0;
     }
   }
@@ -695 +766 @@
     G4int k = 0;
     G4double fac = 1.0;
-    if (nmod > 1 && e > upperEkin[regModels->ModelIndex(0)]) {
+    if (nmod > 1 && e > eLow[1]) {
       do {
         k++;
-        fac *= (1.0 + factor[k]*upperEkin[regModels->ModelIndex(k-1)]/e);
-      } while (k<nmod-1 && e > upperEkin[regModels->ModelIndex(k)] );
+        fac *= (1.0 + factor[k]/e);
+      } while ( k+1 < nmod && e > eLow[k+1] );
     }
 
@@ -721 +792 @@
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4EmModelManager::DumpModelList(G4int verb)
+{
+  if(verb == 0) return;
+  for(G4int i=0; i<nRegions; i++) {
+    G4RegionModels* r = setOfRegionModels[i];
+    const G4Region* reg = r->Region();
+    if(verb > 1 || nRegions > 1) {
+    }
+    G4int n = r->NumberOfModels();  
+    if(verb > 1 || n > 0) {
+      G4cout << "      ===== EM models for the G4Region  " << reg->GetName()
+             << " ======" << G4endl;;
+      for(G4int j=0; j<n; j++) {
+        const G4VEmModel* m = models[r->ModelIndex(j)];
+        G4cout << std::setw(20);
+        G4cout << m->GetName() << " :     Emin= " 
+               << std::setw(10) << G4BestUnit(r->LowEdgeEnergy(j),"Energy")
+               << "        Emax=   " 
+               << G4BestUnit(r->LowEdgeEnergy(j+1),"Energy")
+               << G4endl;
+      }  
+    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

trunk/source/processes/electromagnetic/utils/src/G4EmMultiModel.cc

@@ -25 +25 @@
 //
 // $Id: G4EmMultiModel.cc,v 1.6 2007/05/22 17:31:58 vnivanch Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------

trunk/source/processes/electromagnetic/utils/src/G4EmProcessOptions.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4EmProcessOptions.cc,v 1.22 2007/11/07 18:38:49 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4EmProcessOptions.cc,v 1.26 2009/02/18 14:43:27 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -59 +59 @@
 #include "G4VMultipleScattering.hh"
 #include "G4Region.hh"
+#include "G4RegionStore.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -392 +393 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4EmProcessOptions::ActivateDeexcitation(G4bool val, const G4Region* r)
-{
-  const std::vector<G4VEnergyLossProcess*>& v =
-        theManager->GetEnergyLossProcessVector();
-  std::vector<G4VEnergyLossProcess*>::const_iterator itr;
-  for(itr = v.begin(); itr != v.end(); itr++) {
-    G4VEnergyLossProcess* p = *itr;
-    if(p) p->ActivateDeexcitation(val,r);
-  }
-  const std::vector<G4VEmProcess*>& w =
-        theManager->GetEmProcessVector();
-  std::vector<G4VEmProcess*>::const_iterator itp;
-  for(itp = w.begin(); itp != w.end(); itp++) {
-    G4VEmProcess* q = *itp;
-    if(q) q->ActivateDeexcitation(val,r);
+void G4EmProcessOptions::ActivateDeexcitation(const G4String& pname, 
+                                              G4bool val, 
+                                              const G4String& reg)
+{
+  G4RegionStore* regionStore = G4RegionStore::GetInstance();
+  const G4Region* r = 0;
+  if(reg == "" || reg == "World") {
+    r = regionStore->GetRegion("DefaultRegionForTheWorld", false);
+  } else {
+    r = regionStore->GetRegion(reg, false);
+  }
+  if(!r) {
+    G4cout << "G4EmProcessOptions::ActivateDeexcitation ERROR: G4Region <"
+           << reg << "> not found, the command ignored" << G4endl;
+    return;
+  }
+
+  const std::vector<G4VEnergyLossProcess*>& v =
+        theManager->GetEnergyLossProcessVector();
+  std::vector<G4VEnergyLossProcess*>::const_iterator itr;
+  for(itr = v.begin(); itr != v.end(); itr++) {
+    G4VEnergyLossProcess* p = *itr;
+    if(p) {
+      if(pname == p->GetProcessName()) p->ActivateDeexcitation(val,r);
+    }
+  }
+  const std::vector<G4VEmProcess*>& w =
+        theManager->GetEmProcessVector();
+  std::vector<G4VEmProcess*>::const_iterator itp;
+  for(itp = w.begin(); itp != w.end(); itp++) {
+    G4VEmProcess* q = *itp;
+    if(q) {
+      if(pname == q->GetProcessName()) q->ActivateDeexcitation(val,r);
+    }
   }
 }
@@ -477 +497 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+void G4EmProcessOptions::SetPolarAngleLimit(G4double val)
+{
+  const std::vector<G4VMultipleScattering*>& u =
+        theManager->GetMultipleScatteringVector();
+  std::vector<G4VMultipleScattering*>::const_iterator itm;
+  for(itm = u.begin(); itm != u.end(); itm++) {
+    if(*itm) (*itm)->SetPolarAngleLimit(val);
+  }
+  const std::vector<G4VEmProcess*>& w =
+        theManager->GetEmProcessVector();
+  std::vector<G4VEmProcess*>::const_iterator itp;
+  for(itp = w.begin(); itp != w.end(); itp++) {
+    G4VEmProcess* q = *itp;
+    if(q) q->SetPolarAngleLimit(val);
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4EmProcessOptions::SetLPMFlag(G4bool val)
 {
@@ -484 +523 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
+void G4EmProcessOptions::SetSplineFlag(G4bool val)
+{
+  theManager->SetSplineFlag(val);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
 void G4EmProcessOptions::SetLinearLossLimit(G4double val)
 {

trunk/source/processes/electromagnetic/utils/src/G4EnergyLossMessenger.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4EnergyLossMessenger.cc,v 1.29 2007/06/11 14:56:51 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4EnergyLossMessenger.cc,v 1.37 2009/02/18 14:43:27 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -49 +49 @@
 // 16-03-07 modify /process/eLoss/minsubsec command (V.Ivanchenko)
 // 18-05-07 add /process/msc directory and commands (V.Ivanchenko)
+// 11-03-08 add /process/em directory and commands (V.Ivanchenko)
 //
 // -------------------------------------------------------------------
@@ -84 +85 @@
   mscDirectory = new G4UIdirectory("/process/msc/");
   mscDirectory->SetGuidance("Commands for EM scattering processes.");
+  emDirectory = new G4UIdirectory("/process/em/");
+  emDirectory->SetGuidance("General commands for EM processes.");
 
   RndmStepCmd = new G4UIcmdWithABool("/process/eLoss/rndmStep",this);
@@ -146 +149 @@
 
   IntegCmd = new G4UIcmdWithABool("/process/eLoss/integral",this);
-  IntegCmd->SetGuidance("Switch true/false the integration of cross section over step.");
+  IntegCmd->SetGuidance("Switch true/false the integral option");
   IntegCmd->SetParameterName("integ",true);
   IntegCmd->SetDefaultValue(true);
@@ -152 +155 @@
 
   rangeCmd = new G4UIcmdWithABool("/process/eLoss/CSDARange",this);
-  rangeCmd->SetGuidance("Switch true/false the precise range calculation.");
+  rangeCmd->SetGuidance("Switch true/false the CSDA range calculation");
   rangeCmd->SetParameterName("range",true);
   rangeCmd->SetDefaultValue(true);
@@ -158 +161 @@
 
   lpmCmd = new G4UIcmdWithABool("/process/eLoss/LPM",this);
-  lpmCmd->SetGuidance("Switch true/false the LPM effect calculation.");
+  lpmCmd->SetGuidance("The flag of the LPM effect calculation");
   lpmCmd->SetParameterName("lpm",true);
   lpmCmd->SetDefaultValue(true);
   lpmCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
+  splCmd = new G4UIcmdWithABool("/process/em/spline",this);
+  splCmd->SetGuidance("The flag of usage spline for Physics Vectors");
+  splCmd->SetParameterName("spl",true);
+  splCmd->SetDefaultValue(false);
+  splCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
+
+  aplCmd = new G4UIcmdWithABool("/process/em/applyCuts",this);
+  aplCmd->SetGuidance("The flag to Apply Cuts for gamma processes");
+  aplCmd->SetParameterName("apl",true);
+  aplCmd->SetDefaultValue(false);
+  aplCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
+
+  deexCmd = new G4UIcommand("/process/em/deexcitation",this);
+  deexCmd->SetGuidance("Set deexcitation flag per process and G4Region.");
+  deexCmd->SetGuidance("  procName  : process name");
+  deexCmd->SetGuidance("  flag      : flag");
+  deexCmd->SetGuidance("  regName   : G4Region name");
+
+  G4UIparameter* pName = new G4UIparameter("pName",'s',false);
+  deexCmd->SetParameter(pName);
+
+  G4UIparameter* flag = new G4UIparameter("flag",'s',false);
+  deexCmd->SetParameter(flag);
+
+  G4UIparameter* regName = new G4UIparameter("regName",'s',false);
+  deexCmd->SetParameter(regName);
+
   dedxCmd = new G4UIcmdWithAnInteger("/process/eLoss/binsDEDX",this);
-  dedxCmd->SetGuidance("Set number of bins for DEDX tables.");
+  dedxCmd->SetGuidance("Set number of bins for DEDX tables");
   dedxCmd->SetParameterName("binsDEDX",true);
   dedxCmd->SetDefaultValue(120);
@@ -170 +200 @@
 
   lamCmd = new G4UIcmdWithAnInteger("/process/eLoss/binsLambda",this);
-  lamCmd->SetGuidance("Set number of bins for Lambda tables.");
+  lamCmd->SetGuidance("Set number of bins for Lambda tables");
   lamCmd->SetParameterName("binsL",true);
   lamCmd->SetDefaultValue(120);
@@ -176 +206 @@
 
   verCmd = new G4UIcmdWithAnInteger("/process/eLoss/verbose",this);
-  verCmd->SetGuidance("Set verbose level for EM physics.");
+  verCmd->SetGuidance("Set verbose level for EM physics");
   verCmd->SetParameterName("verb",true);
   verCmd->SetDefaultValue(1);
   verCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
+  ver1Cmd = new G4UIcmdWithAnInteger("/process/em/verbose",this);
+  ver1Cmd->SetGuidance("Set verbose level for EM physics");
+  ver1Cmd->SetParameterName("verb1",true);
+  ver1Cmd->SetDefaultValue(1);
+  ver1Cmd->AvailableForStates(G4State_PreInit,G4State_Idle);
+
   lllCmd = new G4UIcmdWithADouble("/process/eLoss/linLossLimit",this);
-  lllCmd->SetGuidance("Set linearLossLimit parameter.");
+  lllCmd->SetGuidance("Set linearLossLimit parameter");
   lllCmd->SetParameterName("linlim",true);
   lllCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
-  labCmd = new G4UIcmdWithADouble("/process/eLoss/lambdaFactor",this);
-  labCmd->SetGuidance("Set lambdaFactor parameter.");
+  labCmd = new G4UIcmdWithADouble("/process/eLoss/LambdaFactor",this);
+  labCmd->SetGuidance("Set lambdaFactor parameter for integral option");
   labCmd->SetParameterName("Fl",true);
   labCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
   mscCmd = new G4UIcmdWithAString("/process/msc/StepLimit",this);
-  mscCmd->SetGuidance("Set msc step limitation type.");
+  mscCmd->SetGuidance("Set msc step limitation type");
   mscCmd->SetParameterName("StepLim",true);
   mscCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
   latCmd = new G4UIcmdWithABool("/process/msc/LateralDisplacement",this);
-  latCmd->SetGuidance("Switch true/false sampling of latra dislacent.");
+  latCmd->SetGuidance("Set flag of sampling of lateral displacement");
   latCmd->SetParameterName("lat",true);
   latCmd->SetDefaultValue(true);
@@ -203 +239 @@
 
   frCmd = new G4UIcmdWithADouble("/process/msc/RangeFactor",this);
-  frCmd->SetGuidance("Set RangeFactor parameter for msc process.");
+  frCmd->SetGuidance("Set RangeFactor parameter for msc processes");
   frCmd->SetParameterName("Fr",true);
   frCmd->SetRange("Fr>0");
@@ -210 +246 @@
 
   fgCmd = new G4UIcmdWithADouble("/process/msc/GeomFactor",this);
-  fgCmd->SetGuidance("Set GeomFactor parameter for msc process.");
+  fgCmd->SetGuidance("Set GeomFactor parameter for msc processes");
   fgCmd->SetParameterName("Fg",true);
   fgCmd->SetRange("Fg>0");
@@ -217 +253 @@
 
   skinCmd = new G4UIcmdWithADouble("/process/msc/Skin",this);
-  skinCmd->SetGuidance("Set skin parameter for multiple scattering.");
+  skinCmd->SetGuidance("Set skin parameter for msc processes");
   skinCmd->SetParameterName("skin",true);
   skinCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 
+  angCmd = new G4UIcmdWithADoubleAndUnit("/process/msc/ThetaLimit",this);
+  angCmd->SetGuidance("Set the limit on the polar angle");
+  angCmd->SetParameterName("theta",true);
+  angCmd->SetUnitCategory("Angle");
+  angCmd->AvailableForStates(G4State_PreInit,G4State_Idle);
 }
 
@@ -233 +274 @@
   delete MinSubSecCmd;
   delete StepFuncCmd;
+  delete deexCmd;
   delete eLossDirectory;
   delete mscDirectory;
+  delete emDirectory;
   delete MinEnCmd;
   delete MaxEnCmd;
@@ -240 +283 @@
   delete rangeCmd;
   delete lpmCmd;
+  delete splCmd;
+  delete aplCmd;
   delete latCmd;
   delete verCmd;
+  delete ver1Cmd;
   delete mscCmd;
   delete dedxCmd;
@@ -250 +296 @@
   delete labCmd;
   delete skinCmd;
+  delete angCmd;
 }
 
@@ -289 +336 @@
   }
 
+  if (command == deexCmd) {
+    G4String s1 (""), s2(""), s3("");
+    G4bool b = false;
+    std::istringstream is(newValue);
+    is >> s1 >> s2 >> s3;
+    if(s2 == "true") b = true;
+    opt->ActivateDeexcitation(s1,b,s3);
+  }
+
   if (command == mscCmd) {
     if(newValue == "Minimal") 
@@ -317 +373 @@
   }  
 
-  if (command == IntegCmd) 
+  if (command == IntegCmd) {
     opt->SetIntegral(IntegCmd->GetNewBoolValue(newValue));
-  
+  }
   if (command == rangeCmd) {
     opt->SetBuildCSDARange(rangeCmd->GetNewBoolValue(newValue));
@@ -330 +386 @@
   }
 
+  if (command == splCmd) {
+    opt->SetSplineFlag(splCmd->GetNewBoolValue(newValue));
+    G4UImanager::GetUIpointer()->ApplyCommand("/run/physicsModified");
+  }
+
+  if (command == aplCmd) {
+    opt->SetApplyCuts(aplCmd->GetNewBoolValue(newValue));
+  }
+
   if (command == latCmd) {
     opt->SetMscLateralDisplacement(latCmd->GetNewBoolValue(newValue));
@@ -335 +400 @@
   }
 
-  if (command == verCmd) 
+  if (command == verCmd) {
     opt->SetVerbose(verCmd->GetNewIntValue(newValue));
-
-  if (command == lllCmd) 
+  }
+  if (command == ver1Cmd) {
+    opt->SetVerbose(ver1Cmd->GetNewIntValue(newValue));
+  }
+  if (command == lllCmd) {
     opt->SetLinearLossLimit(lllCmd->GetNewDoubleValue(newValue));
-
-  if (command == labCmd) 
+  }
+  if (command == labCmd) {
     opt->SetLambdaFactor(labCmd->GetNewDoubleValue(newValue));
-  
+  }
   if (command == skinCmd) {
     opt->SetSkin(skinCmd->GetNewDoubleValue(newValue));
@@ -364 +432 @@
     G4UImanager::GetUIpointer()->ApplyCommand("/run/physicsModified");
   }
+  if (command == angCmd) { 
+    opt->SetPolarAngleLimit(angCmd->GetNewDoubleValue(newValue));
+    G4UImanager::GetUIpointer()->ApplyCommand("/run/physicsModified");
+  }  
 }
 

trunk/source/processes/electromagnetic/utils/src/G4EnergyLossTables.cc

@@ -25 +25 @@
 //
 //
-// $Id: G4EnergyLossTables.cc,v 1.33 2006/06/29 19:55:09 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4EnergyLossTables.cc,v 1.34 2008/07/08 10:57:22 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -998 +998 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-void G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition* aParticle, const G4String& q)
-{
-  G4String s = "G4EnergyLossTables:: " + q + " table not found for "
-             + aParticle->GetParticleName() + "!";
-  G4Exception(s);
-  exit(1);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void G4EnergyLossTables::ParticleHaveNoLoss(const G4ParticleDefinition* aParticle, 
+                                            const G4String& q)
+{
+  G4String s = " " + q + " table not found for "
+             + aParticle->GetParticleName() + " !";
+  G4Exception("G4EnergyLossTables::ParticleHaveNoLoss", "EM01",
+              FatalException, s);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

trunk/source/processes/electromagnetic/utils/src/G4LossTableBuilder.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4LossTableBuilder.cc,v 1.24 2007/02/16 11:59:35 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4LossTableBuilder.cc,v 1.28 2009/02/18 16:24:47 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -64 +64 @@
 
 G4LossTableBuilder::G4LossTableBuilder() 
-{}
+{
+  splineFlag = true;
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -73 +75 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
-                                  const std::vector<G4PhysicsTable*>& list)
+void 
+G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
+                                   const std::vector<G4PhysicsTable*>& list)
 {
   size_t n_processes = list.size();
@@ -91 +94 @@
 
     pv = new G4PhysicsLogVector(elow, ehigh, nbins);
+    pv->SetSpline(splineFlag);
     for (size_t j=0; j<nbins; j++) {
       G4double dedx = 0.0;
@@ -128 +132 @@
       G4double dedx1  = pv->GetValue(elow, b);
 
+      //G4cout << "nbins= " << nbins << " dedx1= " << dedx1 << G4endl;
+
+      // protection for specific cases dedx=0
       if(dedx1 == 0.0) {
-        for (size_t k=1; k<nbins; k++) {
+        for (size_t k=1; k<nbins-1; k++) {
           bin0++;
           elow  = pv->GetLowEdgeEnergy(k);
@@ -138 +145 @@
       }
 
+      //G4cout << "nbins= " << nbins << " elow= " << elow << " ehigh= " << ehigh << G4endl;
+      // initialisation of a new vector
       G4PhysicsLogVector* v = new G4PhysicsLogVector(elow, ehigh, nbins);
-
+      // dedx is exect zero
+      if(nbins <= 2) {
+        v->PutValue(0,1000.);
+        v->PutValue(1,2000.);
+        G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
+        return;
+      }
+      v->SetSpline(splineFlag);
+
+      // assumed dedx proportional to beta
       G4double range  = 2.*elow/dedx1;
-      //G4double range  = elow/dedx1;
       v->PutValue(0,range);
       G4double energy1 = elow;
@@ -148 +165 @@
 
         G4double energy2 = pv->GetLowEdgeEnergy(j+bin0);
-        G4double dedx2   = pv->GetValue(energy2, b);
         G4double de      = (energy2 - energy1) * del;
-        G4double energy  = energy1 - de*0.5;
-
-        G4bool   yes     = true;
-        //G4bool   yes     = false;
-        if(dedx1 < DBL_MIN || dedx2 < DBL_MIN) yes = false;   
-
-        G4double fac, f;
-
-        if(yes) fac = std::log(dedx2/dedx1)/std::log(energy2/energy1);
-        else    fac = (dedx2 - dedx1)/(energy2 - energy1);
-
+        G4double energy  = energy2 + de*0.5;
         for (size_t k=0; k<n; k++) {
-          energy += de;
-          if(yes) f = dedx1*std::exp(fac*std::log(energy/energy1));
-          else    f = dedx1 + fac*(energy - energy1);
-          if(f > DBL_MIN) range  += de/f;
-        }
+          energy -= de;
+          dedx1 = pv->GetValue(energy, b);
+          if(dedx1 > 0.0) range += de/dedx1;
+        }
+
         //      G4cout << "Range i= " <<i << " j= " << j << G4endl;
         v->PutValue(j,range);
         energy1 = energy2;
-        dedx1   = dedx2;
       }
       G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
@@ -197 +202 @@
       G4double rlow  = pv->GetValue(elow, b);
       G4double rhigh = pv->GetValue(ehigh, b);
-
-      rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
-
       
       G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(nbins,rlow,rhigh);
+      v->SetSpline(splineFlag);
+
       for (size_t j=0; j<nbins; j++) {
         G4double e  = pv->GetLowEdgeEnergy(j);
@@ -207 +211 @@
         v->PutValues(j,r,e);
       }
+      v->PutValues(nbins,rhigh+rlow,ehigh);
 
       G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);

trunk/source/processes/electromagnetic/utils/src/G4LossTableManager.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4LossTableManager.cc,v 1.84 2007/06/14 07:28:48 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4LossTableManager.cc,v 1.95 2008/11/13 18:23:39 schaelic Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -70 +70 @@
 // 12-02-07 Add SetSkin, SetLinearLossLimit (V.Ivanchenko)
 // 18-06-07 Move definition of msc parameters to G4EmProcessOptions (V.Ivanchenko)
+// 21-02-08 Add G4EmSaturation (V.Ivanchenko)
 //
 // Class Description:
@@ -91 +92 @@
 #include "G4PhysicsTableHelper.hh"
 #include "G4EmCorrections.hh"
+#include "G4EmSaturation.hh"
 #include "G4EmTableType.hh"
 #include "G4LossTableBuilder.hh"
@@ -116 +118 @@
   size_t msc = msc_vector.size();
   for (size_t j=0; j<msc; j++) {
-    if(msc_vector[j] ) delete msc_vector[j];
+    if( msc_vector[j] ) delete msc_vector[j];
   }
   size_t emp = emp_vector.size();
   for (size_t k=0; k<emp; k++) {
-    if(emp_vector[k] ) delete emp_vector[k];
+    if( emp_vector[k] ) delete emp_vector[k];
+  }
+  size_t mod = mod_vector.size();
+  for (size_t a=0; a<mod; a++) {
+    if( mod_vector[a] ) delete mod_vector[a];
+  }
+  size_t fmod = fmod_vector.size();
+  for (size_t b=0; b<fmod; b++) {
+    if( fmod_vector[b] ) delete fmod_vector[b];
   }
   Clear();
@@ -126 +136 @@
   delete tableBuilder;
   delete emCorrections;
+  delete emSaturation;
 }
 
@@ -152 +163 @@
   tableBuilder = new G4LossTableBuilder();
   emCorrections= new G4EmCorrections();
+  emSaturation = new G4EmSaturation();
   integral = true;
   integralActive = false;
@@ -160 +172 @@
   stepFunctionActive = false;
   flagLPM = true;
+  splineFlag = true;
   bremsTh = DBL_MAX;
   verbose = 1;
+  tableBuilder->SetSplineFlag(splineFlag);
 }
 
@@ -226 +240 @@
 {
   msc_vector.push_back(p);
-  if(verbose > 1) 
+  if(verbose > 1) {
     G4cout << "G4LossTableManager::Register G4VMultipleScattering : " 
            << p->GetProcessName() << G4endl;
+  }
 }
 
@@ -246 +261 @@
 {
   emp_vector.push_back(p);
-  if(verbose > 1) 
+  if(verbose > 1) {
     G4cout << "G4LossTableManager::Register G4VEmProcess : " 
            << p->GetProcessName() << G4endl;
+  }
 }
 
@@ -258 +274 @@
   for (size_t i=0; i<emp; i++) {
     if(emp_vector[i] == p) emp_vector[i] = 0;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+void G4LossTableManager::Register(G4VEmModel* p)
+{
+  mod_vector.push_back(p);
+  if(verbose > 1) {
+    G4cout << "G4LossTableManager::Register G4VEmModel : " 
+           << p->GetName() << G4endl;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+void G4LossTableManager::DeRegister(G4VEmModel* p)
+{
+  size_t n = mod_vector.size();
+  for (size_t i=0; i<n; i++) {
+    if(mod_vector[i] == p) mod_vector[i] = 0;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+void G4LossTableManager::Register(G4VEmFluctuationModel* p)
+{
+  fmod_vector.push_back(p);
+  if(verbose > 1) {
+    G4cout << "G4LossTableManager::Register G4VEmFluctuationModel : " 
+           << p->GetName() << G4endl;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+void G4LossTableManager::DeRegister(G4VEmFluctuationModel* p)
+{
+  size_t n = fmod_vector.size();
+  for (size_t i=0; i<n; i++) {
+    if(fmod_vector[i] == p) fmod_vector[i] = 0;
   }
 }
@@ -364 +422 @@
      const G4ParticleDefinition* aParticle)
 {
-  G4String s = "G4LossTableManager:: dE/dx table not found for "
-             + aParticle->GetParticleName() + "!";
-  G4Exception(s);
-  exit(1);
+  G4String s = " dE/dx table not found for "
+             + aParticle->GetParticleName() + " !";
+  G4Exception("G4LossTableManager::ParticleHaveNoLoss", "EM01",
+              FatalException, s);
+
 }
 
@@ -476 +535 @@
     p = loss_vector[i];
     if (p && aParticle == part_vector[i] && !tables_are_built[i]) {
-      if (p->IsIonisationProcess() && isActive[i] || !em || (em && !isActive[iem]) ) {
+      if ((p->IsIonisationProcess() && isActive[i]) || 
+          !em || (em && !isActive[iem]) ) {
         em = p;
         iem= i;
@@ -509 +569 @@
   dedx = em->IonisationTable();
   if (1 < n_dedx) {
-    em->SetDEDXTable(dedx, fIonisation);
+    em->SetDEDXTable(dedx, fIsIonisation);
     dedx = 0;
     dedx  = G4PhysicsTableHelper::PreparePhysicsTable(dedx);
@@ -560 +620 @@
     G4PhysicsTable* dedxSub = em->IonisationTableForSubsec();
     if (1 < listSub.size()) {
-      em->SetDEDXTable(dedxSub, fSubIonisation);
+      em->SetDEDXTable(dedxSub, fIsSubIonisation);
       dedxSub = 0;
       dedxSub = G4PhysicsTableHelper::PreparePhysicsTable(dedxSub);
@@ -829 +889 @@
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+void G4LossTableManager::SetSplineFlag(G4bool val)
+{
+  splineFlag = val;
+  tableBuilder->SetSplineFlag(splineFlag);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+G4bool G4LossTableManager::SplineFlag() const
+{
+  return splineFlag;
+}
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
@@ -843 +918 @@
 }
 
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4EmCorrections* G4LossTableManager::EmCorrections() 
+{
+  return emCorrections;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4EmSaturation* G4LossTableManager::EmSaturation()
+{
+  return emSaturation;
+}
+
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....

trunk/source/processes/electromagnetic/utils/src/G4VEmFluctuationModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEmFluctuationModel.cc,v 1.2 2006/06/29 19:55:15 gunter Exp $
-// GEANT4 tag $Name:  $
+// $Id: G4VEmFluctuationModel.cc,v 1.4 2009/02/19 11:25:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -49 +49 @@
 
 #include "G4VEmFluctuationModel.hh"
-
+#include "G4LossTableManager.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -56 +56 @@
 G4VEmFluctuationModel::G4VEmFluctuationModel(const G4String& nam)
   : name(nam) 
+{
+  G4LossTableManager::Instance()->Register(this);
+}
+
+G4VEmFluctuationModel::~G4VEmFluctuationModel() 
+{
+  G4LossTableManager::Instance()->DeRegister(this);
+}
+
+void G4VEmFluctuationModel::InitialiseMe(const G4ParticleDefinition*)
 {}
 
-G4VEmFluctuationModel::~G4VEmFluctuationModel() 
+void G4VEmFluctuationModel::SetParticleAndCharge(const G4ParticleDefinition*, 
+                                                 G4double)
 {}
 

trunk/source/processes/electromagnetic/utils/src/G4VEmModel.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEmModel.cc,v 1.8 2007/09/25 10:19:07 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4VEmModel.cc,v 1.25 2009/02/19 09:57:36 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -41 +41 @@
 // 25.10.2005 Set default highLimit=100.TeV (V.Ivanchenko)
 // 06.02.2006 add method ComputeMeanFreePath() (mma)
+// 16.02.2009 Move implementations of virtual methods to source (VI)
 //
 //
@@ -51 +52 @@
 
 #include "G4VEmModel.hh"
+#include "G4LossTableManager.hh"
+#include "G4ProductionCutsTable.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -56 +59 @@
 
 G4VEmModel::G4VEmModel(const G4String& nam):
- lowLimit(0.1*keV), highLimit(100.0*TeV), fluc(0), name(nam), pParticleChange(0)
-{}
+  fluc(0), name(nam), lowLimit(0.1*keV), highLimit(100.0*TeV), 
+  polarAngleLimit(0.0),secondaryThreshold(DBL_MAX),theLPMflag(false),
+  pParticleChange(0),nuclearStopping(false),
+  currentCouple(0),currentElement(0),
+  nsec(5),flagDeexcitation(false) 
+{
+  xsec.resize(nsec);
+  nSelectors = 0;
+  G4LossTableManager::Instance()->Register(this);
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
 G4VEmModel::~G4VEmModel()
-{}
+{
+  G4LossTableManager::Instance()->DeRegister(this);
+  G4int n = elmSelectors.size();
+  if(n > 0) {
+    for(G4int i=0; i<n; i++) { 
+      delete elmSelectors[i]; 
+    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p, 
+                                            const G4DataVector& cuts)
+{
+  // initialise before run
+  flagDeexcitation = false;
+
+  G4int nbins = G4int(std::log10(highLimit/lowLimit) + 0.5);
+  if(nbins < 3) nbins = 3;
+  G4bool spline = G4LossTableManager::Instance()->SplineFlag();
+
+  G4ProductionCutsTable* theCoupleTable=
+    G4ProductionCutsTable::GetProductionCutsTable();
+  G4int numOfCouples = theCoupleTable->GetTableSize();
+
+  // prepare vector
+  if(numOfCouples > nSelectors) {
+    elmSelectors.resize(numOfCouples);
+    nSelectors = numOfCouples;
+  }
+
+  // initialise vector
+  for(G4int i=0; i<numOfCouples; i++) {
+    const G4MaterialCutsCouple* couple =
+      theCoupleTable->GetMaterialCutsCouple(i);
+    const G4Material* material = couple->GetMaterial();
+    G4int idx = couple->GetIndex();
+
+    // selector already exist check if should be deleted
+    G4bool create = true;
+    if(elmSelectors[i]) {
+      if(material == elmSelectors[i]->GetMaterial()) create = false;
+      else delete elmSelectors[i];
+    }
+    if(create) {
+      elmSelectors[i] = new G4EmElementSelector(this,material,nbins,
+                                                lowLimit,highLimit,spline);
+    }
+    elmSelectors[i]->Initialise(p, cuts[idx]);
+    //elmSelectors[i]->Dump(p);
+  } 
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
+                                          const G4ParticleDefinition*,
+                                          G4double,G4double)
+{
+  return 0.0;
+}
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
@@ -68 +140 @@
 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
                                            const G4ParticleDefinition* p,
-                                                 G4double ekin,
-                                                 G4double emin,
-                                                 G4double emax)
-{
+                                           G4double ekin,
+                                           G4double emin,
+                                           G4double emax)
+{
+  SetupForMaterial(p, material, ekin);
   G4double cross = 0.0;
   const G4ElementVector* theElementVector = material->GetElementVector();
   const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
-  size_t nelm = material->GetNumberOfElements();
-  for (size_t i=0; i<nelm; i++) {
-    const G4Element* elm = (*theElementVector)[i];
+  G4int nelm = material->GetNumberOfElements(); 
+  if(nelm > nsec) {
+    xsec.resize(nelm);
+    nsec = nelm;
+  }
+  for (G4int i=0; i<nelm; i++) {
     cross += theAtomNumDensityVector[i]*
-      ComputeCrossSectionPerAtom(p,ekin,elm->GetZ(),elm->GetN(),emin,emax);
+      ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
     xsec[i] = cross;
   }
@@ -87 +163 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 
-G4double G4VEmModel::ComputeMeanFreePath(const G4ParticleDefinition* p,
-                                               G4double ekin,
-                                         const G4Material* material,     
-                                               G4double emin,
-                                               G4double emax)
-{
-  G4double mfp = DBL_MAX;
-  G4double cross = CrossSectionPerVolume(material,p,ekin,emin,emax);
-  if (cross > DBL_MIN) mfp = 1./cross;
-  return mfp;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
-
-
+G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
+                                                G4double, G4double, G4double,
+                                                G4double, G4double)
+{
+  return 0.0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::MinEnergyCut(const G4ParticleDefinition*,
+                                  const G4MaterialCutsCouple*)
+{
+  return 0.0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
+                                          const G4Material*, G4double)
+{
+  G4double q = p->GetPDGCharge()/CLHEP::eplus;
+  return q*q;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
+                                       const G4Material*, G4double)
+{
+  return p->GetPDGCharge();
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
+                                      const G4DynamicParticle*,
+                                      G4double&,G4double&,G4double)
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4VEmModel::SampleDeexcitationAlongStep(const G4Material*,
+                                             const G4Track&,
+                                             G4double& )
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
+                                        G4double kineticEnergy)
+{
+  return kineticEnergy;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4VEmModel::SampleScattering(const G4DynamicParticle*, G4double)
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::ComputeTruePathLengthLimit(const G4Track&, 
+                                                G4PhysicsTable*, 
+                                                G4double)
+{
+  return DBL_MAX;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::ComputeGeomPathLength(G4double truePathLength)
+{
+  return truePathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+G4double G4VEmModel::ComputeTrueStepLength(G4double geomPathLength)
+{
+  return geomPathLength;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4VEmModel::DefineForRegion(const G4Region*) 
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+
+void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
+                                  const G4Material*, G4double)
+{}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

trunk/source/processes/electromagnetic/utils/src/G4VEmProcess.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEmProcess.cc,v 1.48 2007/10/29 08:38:58 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4VEmProcess.cc,v 1.62 2009/02/19 09:57:36 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -82 +82 @@
 
 G4VEmProcess::G4VEmProcess(const G4String& name, G4ProcessType type):
-                      G4VDiscreteProcess(name, type),
-  selectedModel(0),                   
+  G4VDiscreteProcess(name, type),
+  secondaryParticle(0),
+  buildLambdaTable(true),
   theLambdaTable(0),
   theEnergyOfCrossSectionMax(0),
   theCrossSectionMax(0),
-  particle(0),
-  secondaryParticle(0),
-  nLambdaBins(90),
-  lambdaFactor(0.8),
-  currentCouple(0),
   integral(false),
-  buildLambdaTable(true),
   applyCuts(false),
   startFromNull(true),
-  nRegions(0)
+  useDeexcitation(false),
+  nDERegions(0),
+  idxDERegions(0),
+  currentModel(0),
+  particle(0),
+  currentCouple(0)
 {
   SetVerboseLevel(1);
+
+  // Size of tables assuming spline
   minKinEnergy = 0.1*keV;
-  maxKinEnergy = 100.0*GeV;
+  maxKinEnergy = 100.0*TeV;
+  nLambdaBins  = 84;
+
+  // default lambda factor
+  lambdaFactor  = 0.8;
+
+  // default limit on polar angle
+  polarAngleLimit = 0.0;
+
+  // particle types
   theGamma     = G4Gamma::Gamma();
   theElectron  = G4Electron::Electron();
@@ -126 +137 @@
   delete modelManager;
   (G4LossTableManager::Instance())->DeRegister(this);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEmProcess::Clear()
+{
+  delete [] theEnergyOfCrossSectionMax;
+  delete [] theCrossSectionMax;
+  delete [] idxDERegions;
+  theEnergyOfCrossSectionMax = 0;
+  theCrossSectionMax = 0;
+  idxDERegions = 0;
+  currentCouple = 0;
+  preStepLambda = 0.0;
+  mfpKinEnergy  = DBL_MAX;
+  deRegions.clear();
+  nDERegions = 0;
 }
 
@@ -153 +181 @@
       theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
   }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4VEmProcess::Clear()
-{
-  if(theEnergyOfCrossSectionMax) delete [] theEnergyOfCrossSectionMax;
-  if(theCrossSectionMax) delete [] theCrossSectionMax;
-  theEnergyOfCrossSectionMax = 0;
-  theCrossSectionMax = 0;
-  currentCouple = 0;
-  preStepLambda = 0.0;
-  mfpKinEnergy  = DBL_MAX;
+  // Sub Cutoff and Deexcitation
+  if (nDERegions>0) {
+
+    const G4ProductionCutsTable* theCoupleTable=
+          G4ProductionCutsTable::GetProductionCutsTable();
+    size_t numOfCouples = theCoupleTable->GetTableSize();
+
+    idxDERegions = new G4bool[numOfCouples];
+
+    for (size_t j=0; j<numOfCouples; j++) {
+
+      const G4MaterialCutsCouple* couple =
+        theCoupleTable->GetMaterialCutsCouple(j);
+      const G4ProductionCuts* pcuts = couple->GetProductionCuts();
+      G4bool reg = false;
+      for(G4int i=0; i<nDERegions; i++) {
+        if(deRegions[i]) {
+          if(pcuts == deRegions[i]->GetProductionCuts()) reg = true;
+        }
+      }
+      idxDERegions[j] = reg;
+    }
+  }
+  if (1 < verboseLevel && nDERegions>0) {
+    G4cout << " Deexcitation is activated for regions: " << G4endl;
+    for (G4int i=0; i<nDERegions; i++) {
+      const G4Region* r = deRegions[i];
+      G4cout << "           " << r->GetName() << G4endl;
+    }
+  }
 }
 
@@ -228 +273 @@
       G4cout << *theLambdaTable << G4endl;
     }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEmProcess::PrintInfoDefinition()
+{
+  if(verboseLevel > 0) {
+    G4cout << G4endl << GetProcessName() << ":   for  "
+           << particle->GetParticleName();
+    if(integral) G4cout << ", integral: 1 ";
+    if(applyCuts) G4cout << ", applyCuts: 1 ";
+    G4cout << "    SubType= " << GetProcessSubType() << G4endl;
+    if(buildLambdaTable) {
+      G4cout << "      Lambda tables from "
+             << G4BestUnit(minKinEnergy,"Energy") 
+             << " to "
+             << G4BestUnit(maxKinEnergy,"Energy")
+             << " in " << nLambdaBins << " bins, spline: " 
+             << (G4LossTableManager::Instance())->SplineFlag()
+             << G4endl;
+    }
+    PrintInfo();
+    modelManager->DumpModelList(verboseLevel);
+  }
+
+  if(verboseLevel > 2 && buildLambdaTable) {
+    G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
+    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
   }
 }
@@ -295 +369 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
-                                       G4double,
-                                       G4ForceCondition* condition)
-{
-  *condition = NotForced;
-  return G4VEmProcess::MeanFreePath(track);
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 G4VParticleChange* G4VEmProcess::PostStepDoIt(const G4Track& track,
                                               const G4Step&)
@@ -333 +397 @@
   }
 
-  G4VEmModel* currentModel = SelectModel(finalT);
-
+  SelectModel(finalT);
+  if(useDeexcitation) {
+    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
+  }
   /*  
   if(0 < verboseLevel) {
@@ -373 +439 @@
 
         } else if (p == thePositron) {
-          if (e < (*theCutsPositron)[currentMaterialIndex]) {
+          if (electron_mass_c2 < (*theCutsGamma)[currentMaterialIndex] &&
+              e < (*theCutsPositron)[currentMaterialIndex]) {
             good = false;
             e += 2.0*electron_mass_c2;
@@ -394 +461 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4VEmProcess::PrintInfoDefinition()
-{
-  if(verboseLevel > 0) {
-    G4cout << G4endl << GetProcessName() << ": " ;
-    PrintInfo();
-    if(integral) {
-      G4cout << "      Integral mode is used  "<< G4endl;
-    }
-  }
-
-  if (!buildLambdaTable)  return;
-  
-  if(verboseLevel > 0) {
-    G4cout << "      tables are built for  "
-           << particle->GetParticleName()
-           << G4endl
-           << "      Lambda tables from "
-           << G4BestUnit(minKinEnergy,"Energy") 
-           << " to "
-           << G4BestUnit(maxKinEnergy,"Energy")
-           << " in " << nLambdaBins << " bins."
-           << G4endl;
-  }
-
-  if(verboseLevel > 1) {
-    G4cout << "Tables are built for " << particle->GetParticleName()
-           << G4endl;
-
-  if(verboseLevel > 2) {
-    G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
-    if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
-    }
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4VEmProcess::MicroscopicCrossSection(G4double kineticEnergy,
-                                         const G4MaterialCutsCouple* couple)
-{
-  // Cross section per atom is calculated
-  DefineMaterial(couple);
-  G4double cross = 0.0;
-  G4bool b;
-  if(theLambdaTable) {
-    cross = (((*theLambdaTable)[currentMaterialIndex])->
-                           GetValue(kineticEnergy, b));
-
-    cross /= currentMaterial->GetTotNbOfAtomsPerVolume();
-  } else {
-    G4VEmModel* model = SelectModel(kineticEnergy);
-    cross = 
-      model->CrossSectionPerVolume(currentMaterial,particle,kineticEnergy);
-  }
-
-  return cross;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
 G4bool G4VEmProcess::StorePhysicsTable(const G4ParticleDefinition* part,
                                        const G4String& directory,
@@ -508 +515 @@
              << G4endl;
     }
+    if((G4LossTableManager::Instance())->SplineFlag()) {
+      size_t n = theLambdaTable->length();
+      for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
+    }
   } else {
     if (1 < verboseLevel) {
@@ -517 +528 @@
 
   return yes;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEmProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
+{
+  G4RegionStore* regionStore = G4RegionStore::GetInstance();
+  const G4Region* reg = r;
+  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
+
+  // the region is in the list
+  if (nDERegions) {
+    for (G4int i=0; i<nDERegions; i++) {
+      if (reg == deRegions[i]) {
+        if(!val) deRegions[i] = 0;
+        return;
+      }
+    }
+  }
+
+  // new region 
+  if(val) {
+    useDeexcitation = true;
+    deRegions.push_back(reg);
+    nDERegions++;
+  } else {
+    useDeexcitation = false;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEmProcess::CrossSectionPerVolume(G4double kineticEnergy,
+                                             const G4MaterialCutsCouple* couple)
+{
+  // Cross section per atom is calculated
+  DefineMaterial(couple);
+  G4double cross = 0.0;
+  G4bool b;
+  if(theLambdaTable) {
+    cross = (((*theLambdaTable)[currentMaterialIndex])->
+                           GetValue(kineticEnergy, b));
+  } else {
+    SelectModel(kineticEnergy);
+    cross = currentModel->CrossSectionPerVolume(currentMaterial,
+                                                particle,kineticEnergy);
+  }
+
+  return cross;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEmProcess::GetMeanFreePath(const G4Track& track,
+                                       G4double,
+                                       G4ForceCondition* condition)
+{
+  *condition = NotForced;
+  return G4VEmProcess::MeanFreePath(track);
 }
 
@@ -568 +638 @@
   G4PhysicsVector* v = 
     new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nLambdaBins);
+  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
   return v;
 }

trunk/source/processes/electromagnetic/utils/src/G4VEnergyLoss.cc

@@ -26 +26 @@
 //
 // $Id: G4VEnergyLoss.cc,v 1.46 2006/06/29 19:55:21 gunter Exp $
-// GEANT4 tag $Name:  $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 

trunk/source/processes/electromagnetic/utils/src/G4VEnergyLossProcess.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VEnergyLossProcess.cc,v 1.123 2008/01/11 19:55:29 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4VEnergyLossProcess.cc,v 1.146 2009/02/19 11:25:50 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -146 +146 @@
 
 G4VEnergyLossProcess::G4VEnergyLossProcess(const G4String& name, 
-  G4ProcessType type): G4VContinuousDiscreteProcess(name, type),
+                                           G4ProcessType type): 
+  G4VContinuousDiscreteProcess(name, type),
+  secondaryParticle(0),
   nSCoffRegions(0),
+  nDERegions(0),
   idxSCoffRegions(0),
+  idxDERegions(0),
   nProcesses(0),
   theDEDXTable(0),
@@ -165 +169 @@
   theEnergyOfCrossSectionMax(0),
   theCrossSectionMax(0),
-  particle(0),
   baseParticle(0),
-  secondaryParticle(0),
-  currentCouple(0),
-  nBins(120),
-  nBinsCSDA(70),
-  nWarnings(0),
-  linLossLimit(0.05),
   minSubRange(0.1),
-  lambdaFactor(0.8),
-  mfpKinEnergy(0.0),
   lossFluctuationFlag(true),
   rndmStepFlag(false),
   tablesAreBuilt(false),
   integral(true),
+  isIon(false),
   isIonisation(true),
-  useSubCutoff(false)
+  useSubCutoff(false),
+  useDeexcitation(false),
+  particle(0),
+  currentCouple(0),
+  nWarnings(0),
+  mfpKinEnergy(0.0)
 {
   SetVerboseLevel(1);
 
-  // Size of tables
+  // low energy limit
   lowestKinEnergy  = 1.*eV;
+
+  // Size of tables assuming spline
   minKinEnergy     = 0.1*keV;
   maxKinEnergy     = 100.0*TeV;
+  nBins            = 84;
   maxKinEnergyCSDA = 1.0*GeV;
+  nBinsCSDA        = 35;
+
+  // default linear loss limit for spline
+  linLossLimit  = 0.01;
 
   // default dRoverRange and finalRange
   SetStepFunction(0.2, 1.0*mm);
 
-  theElectron = G4Electron::Electron();
-  thePositron = G4Positron::Positron();
+  // default lambda factor
+  lambdaFactor  = 0.8;
+
+  // particle types
+  theElectron   = G4Electron::Electron();
+  thePositron   = G4Positron::Positron();
+  theGenericIon = 0;
 
   // run time objects
@@ -208 +221 @@
   fluctModel = 0;
 
-  scoffRegions.clear();
-  scProcesses.clear();
   scTracks.reserve(5);
   secParticles.reserve(5);
@@ -215 +226 @@
   // Data for stragling of ranges from ICRU'37 report
   const G4int nrbins = 7;
-  vstrag = new G4PhysicsLogVector(keV, GeV, nrbins);
+  vstrag = new G4PhysicsLogVector(keV, GeV, nrbins-1);
+  vstrag->SetSpline(true);
   G4double s[nrbins] = {-0.2, -0.85, -1.3, -1.578, -1.76, -1.85, -1.9};
   for(G4int i=0; i<nrbins; i++) {vstrag->PutValue(i, s[i]);}
@@ -228 +240 @@
            << G4endl;
   delete vstrag;
-  Clear();
+  Clean();
 
   if ( !baseParticle ) {
@@ -282 +294 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4VEnergyLossProcess::Clear()
-{
-  if(1 < verboseLevel) 
+void G4VEnergyLossProcess::Clean()
+{
+  if(1 < verboseLevel) { 
     G4cout << "G4VEnergyLossProcess::Clear() for " << GetProcessName() 
-    << G4endl;
-
+           << G4endl;
+  }
   delete [] theDEDXAtMaxEnergy;
   delete [] theRangeAtMaxEnergy;
@@ -293 +305 @@
   delete [] theCrossSectionMax;
   delete [] idxSCoffRegions;
+  delete [] idxDERegions;
 
   theDEDXAtMaxEnergy = 0;
@@ -300 +313 @@
   tablesAreBuilt = false;
 
-  scTracks.clear();
+  //scTracks.clear();
   scProcesses.clear();
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4VEnergyLossProcess::PreparePhysicsTable(
-     const G4ParticleDefinition& part)
-{
-
-  // Are particle defined?
-  if( !particle ) {
-    if(part.GetParticleType() == "nucleus" && 
-       part.GetParticleSubType() == "generic") 
-         particle = G4GenericIon::GenericIon();
-    else particle = &part;
-  }
-
+  nProcesses = 0;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEnergyLossProcess::MinPrimaryEnergy(const G4ParticleDefinition*, 
+                                                const G4Material*, 
+                                                G4double cut)
+{
+  return cut;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void 
+G4VEnergyLossProcess::PreparePhysicsTable(const G4ParticleDefinition& part)
+{
   if(1 < verboseLevel) {
     G4cout << "G4VEnergyLossProcess::PreparePhysicsTable for "
            << GetProcessName()
            << " for " << part.GetParticleName()
-           << " local: " << particle->GetParticleName()
            << G4endl;
   }
-
-  G4LossTableManager* lManager = G4LossTableManager::Instance();
-
-  if(&part != particle) {
-    if(part.GetParticleType() == "nucleus") lManager->RegisterIon(&part, this);
-    else                          lManager->RegisterExtraParticle(&part, this);
-    return;
-  }
-
-  Clear();
 
   currentCouple = 0;
@@ -341 +344 @@
   fRange        = DBL_MAX;
   preStepKinEnergy = 0.0;
+  chargeSqRatio = 1.0;
+  massRatio = 1.0;
+  reduceFactor = 1.0;
+
+  G4LossTableManager* lManager = G4LossTableManager::Instance();
+
+  // Are particle defined?
+  if( !particle ) {
+    particle = &part;
+    if(part.GetParticleType() == "nucleus") {
+      if(!theGenericIon) theGenericIon = G4GenericIon::GenericIon();
+      if(particle == theGenericIon) { isIon = true; }
+      else if(part.GetPDGCharge() > eplus) {
+        isIon = true; 
+
+        // generic ions created on-fly
+        if(part.GetPDGCharge() > 2.5*eplus) {
+          particle = theGenericIon;
+        }
+      }
+    }
+  }
+
+  if( particle != &part) {
+    if(part.GetParticleType() == "nucleus") {
+      isIon = true;
+      lManager->RegisterIon(&part, this);
+    } else { 
+      lManager->RegisterExtraParticle(&part, this);
+    }
+    return;
+  }
+
+  Clean();
 
   // Base particle and set of models can be defined here
@@ -372 +409 @@
   G4double initialCharge = particle->GetPDGCharge();
   G4double initialMass   = particle->GetPDGMass();
-  chargeSquare = initialCharge*initialCharge/(eplus*eplus);
-  chargeSqRatio = 1.0;
-  massRatio = 1.0;
-  reduceFactor = 1.0;
 
   if (baseParticle) {
@@ -387 +420 @@
                                      minSubRange, verboseLevel);
 
-  // Sub Cutoff Regime
-  scProcesses.clear();
-  nProcesses = 0;
-  
-  if (nSCoffRegions>0) {
+  // Sub Cutoff and Deexcitation
+  if (nSCoffRegions>0 || nDERegions>0) {
     theSubCuts = modelManager->SubCutoff();
 
@@ -397 +427 @@
           G4ProductionCutsTable::GetProductionCutsTable();
     size_t numOfCouples = theCoupleTable->GetTableSize();
-    idxSCoffRegions = new G4int[numOfCouples];
+
+    if(nSCoffRegions>0) idxSCoffRegions = new G4bool[numOfCouples];
+    if(nDERegions>0) idxDERegions = new G4bool[numOfCouples];
   
     for (size_t j=0; j<numOfCouples; j++) {
@@ -404 +436 @@
         theCoupleTable->GetMaterialCutsCouple(j);
       const G4ProductionCuts* pcuts = couple->GetProductionCuts();
-      G4int reg = 0;
-      for(G4int i=0; i<nSCoffRegions; i++) {
-        if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = 1;
-      }
-      idxSCoffRegions[j] = reg;
+      
+      if(nSCoffRegions>0) {
+        G4bool reg = false;
+        for(G4int i=0; i<nSCoffRegions; i++) {
+          if( pcuts == scoffRegions[i]->GetProductionCuts()) reg = true;
+        }
+        idxSCoffRegions[j] = reg;
+      }
+      if(nDERegions>0) {
+        G4bool reg = false;
+        for(G4int i=0; i<nDERegions; i++) {
+          if( pcuts == deRegions[i]->GetProductionCuts()) reg = true;
+        }
+        idxDERegions[j] = reg;
+      }
     }
   }
@@ -416 +458 @@
   if (1 < verboseLevel) {
     G4cout << "G4VEnergyLossProcess::Initialise() is done "
+           << " for local " << particle->GetParticleName()
+           << " isIon= " << isIon
            << " chargeSqRatio= " << chargeSqRatio
            << " massRatio= " << massRatio
@@ -426 +470 @@
       }
     }
+    if (nDERegions) {
+      G4cout << " Deexcitation is ON for regions: " << G4endl;
+      for (G4int i=0; i<nDERegions; i++) {
+        const G4Region* r = deRegions[i];
+        G4cout << "           " << r->GetName() << G4endl;
+      }
+    }
   }
 }
@@ -442 +493 @@
   }
 
-  if(!tablesAreBuilt && &part == particle)
-    G4LossTableManager::Instance()->BuildPhysicsTable(particle, this);
-
-  if(0 < verboseLevel && (&part == particle) && !baseParticle) {
-    PrintInfoDefinition();
-    safetyHelper->InitialiseHelper();
+  if(&part == particle) {
+    if(!tablesAreBuilt) {
+      G4LossTableManager::Instance()->BuildPhysicsTable(particle, this);
+    }
+    if(!baseParticle) {
+      if(0 < verboseLevel) PrintInfoDefinition();
+    
+      // needs to be done only once
+      safetyHelper->InitialiseHelper();
+    }
   }
 
@@ -456 +511 @@
     if(isIonisation) G4cout << "  isIonisation  flag = 1";
     G4cout << G4endl;
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
-{
-  G4RegionStore* regionStore = G4RegionStore::GetInstance();
-  if(val) {
-    useSubCutoff = true;
-    if (!r) r = regionStore->GetRegion("DefaultRegionForTheWorld", false);
-    if (nSCoffRegions) {
-      for (G4int i=0; i<nSCoffRegions; i++) {
-        if (r == scoffRegions[i]) return;
-      }
-    }
-    scoffRegions.push_back(r);
-    nSCoffRegions++;
-  } else {
-    useSubCutoff = false;
   }
 }
@@ -528 +563 @@
   for(size_t i=0; i<numOfCouples; i++) {
 
-    if(1 < verboseLevel) 
+    if(1 < verboseLevel) {
       G4cout << "G4VEnergyLossProcess::BuildDEDXVector flag=  " 
              << table->GetFlag(i) << G4endl;
-
+    }
     if (table->GetFlag(i)) {
 
@@ -538 +573 @@
         theCoupleTable->GetMaterialCutsCouple(i);
       G4PhysicsVector* aVector = new G4PhysicsLogVector(emin, emax, bin);
+      aVector->SetSpline((G4LossTableManager::Instance())->SplineFlag());
+
       modelManager->FillDEDXVector(aVector, couple, tType);
 
@@ -580 +617 @@
            << G4endl;
   }
-  if(!table) return table;
+  if(!table) {return table;}
 
   // Access to materials
@@ -615 +652 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4double G4VEnergyLossProcess::GetContinuousStepLimit(
-                const G4Track&,
-                G4double, G4double, G4double&)
-{
-  return DBL_MAX;
+void G4VEnergyLossProcess::PrintInfoDefinition()
+{
+  if(0 < verboseLevel) {
+    G4cout << G4endl << GetProcessName() << ":   for  "
+           << particle->GetParticleName()
+           << "    SubType= " << GetProcessSubType() 
+           << G4endl
+           << "      dE/dx and range tables from "
+           << G4BestUnit(minKinEnergy,"Energy")
+           << " to " << G4BestUnit(maxKinEnergy,"Energy")
+           << " in " << nBins << " bins" << G4endl
+           << "      Lambda tables from threshold to "
+           << G4BestUnit(maxKinEnergy,"Energy")
+           << " in " << nBins << " bins, spline: " 
+           << (G4LossTableManager::Instance())->SplineFlag()
+           << G4endl;
+    if(theRangeTableForLoss && isIonisation) {
+      G4cout << "      finalRange(mm)= " << finalRange/mm
+             << ", dRoverRange= " << dRoverRange
+             << ", integral: " << integral
+             << ", fluct: " << lossFluctuationFlag
+             << ", linLossLimit= " << linLossLimit
+             << G4endl;
+    }
+    PrintInfo();
+    modelManager->DumpModelList(verboseLevel);
+    if(theCSDARangeTable && isIonisation) {
+      G4cout << "      CSDA range table up"
+             << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
+             << " in " << nBinsCSDA << " bins" << G4endl;
+    }
+    if(nSCoffRegions>0 && isIonisation) {
+      G4cout << "      Subcutoff sampling in " << nSCoffRegions 
+             << " regions" << G4endl;
+    }
+    if(2 < verboseLevel) {
+      G4cout << "      DEDXTable address= " << theDEDXTable << G4endl;
+      if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
+      G4cout << "non restricted DEDXTable address= " 
+             << theDEDXunRestrictedTable << G4endl;
+      if(theDEDXunRestrictedTable && isIonisation) {
+           G4cout << (*theDEDXunRestrictedTable) << G4endl;
+      }
+      if(theDEDXSubTable && isIonisation) {
+        G4cout << (*theDEDXSubTable) << G4endl;
+      }
+      G4cout << "      CSDARangeTable address= " << theCSDARangeTable 
+             << G4endl;
+      if(theCSDARangeTable && isIonisation) {
+        G4cout << (*theCSDARangeTable) << G4endl;
+      }
+      G4cout << "      RangeTableForLoss address= " << theRangeTableForLoss 
+             << G4endl;
+      if(theRangeTableForLoss && isIonisation) {
+             G4cout << (*theRangeTableForLoss) << G4endl;
+      }
+      G4cout << "      InverseRangeTable address= " << theInverseRangeTable 
+             << G4endl;
+      if(theInverseRangeTable && isIonisation) {
+             G4cout << (*theInverseRangeTable) << G4endl;
+      }
+      G4cout << "      LambdaTable address= " << theLambdaTable << G4endl;
+      if(theLambdaTable && isIonisation) {
+        G4cout << (*theLambdaTable) << G4endl;
+      }
+      G4cout << "      SubLambdaTable address= " << theSubLambdaTable << G4endl;
+      if(theSubLambdaTable && isIonisation) {
+        G4cout << (*theSubLambdaTable) << G4endl;
+      }
+    }
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEnergyLossProcess::ActivateSubCutoff(G4bool val, const G4Region* r)
+{
+  G4RegionStore* regionStore = G4RegionStore::GetInstance();
+  const G4Region* reg = r;
+  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
+
+  // the region is in the list
+  if (nSCoffRegions) {
+    for (G4int i=0; i<nSCoffRegions; i++) {
+      if (reg == scoffRegions[i]) {
+        if(!val) deRegions[i] = 0;
+        return;
+      }
+    }
+  }
+
+  // new region 
+  if(val) {
+    useSubCutoff = true;
+    scoffRegions.push_back(reg);
+    nSCoffRegions++;
+  } else {
+    useSubCutoff = false;
+  }
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+void G4VEnergyLossProcess::ActivateDeexcitation(G4bool val, const G4Region* r)
+{
+  G4RegionStore* regionStore = G4RegionStore::GetInstance();
+  const G4Region* reg = r;
+  if (!reg) {reg = regionStore->GetRegion("DefaultRegionForTheWorld", false);}
+
+  // the region is in the list
+  if (nDERegions) {
+    for (G4int i=0; i<nDERegions; i++) {
+      if (reg == deRegions[i]) {
+        if(!val) deRegions[i] = 0;
+        return;
+      }
+    }
+  }
+
+  // new region 
+  if(val) {
+    useDeexcitation = true;
+    deRegions.push_back(reg);
+    nDERegions++;
+  } else {
+    useDeexcitation = false;
+  }
 }
 
@@ -641 +800 @@
     if(x > finalRange && y < currentMinStep) { 
       x = y + finalRange*(1.0 - dRoverRange)*(2.0 - finalRange/fRange);
-    } else if (rndmStepFlag) x = SampleRange();
-    //    G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
-    //    <<" range= "<<fRange <<" cMinSt="<<currentMinStep
-    //    <<" safety= " << safety<< " limit= " << x <<G4endl;
-  }
-  //  G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
+    } else if (rndmStepFlag) {x = SampleRange();}
+    //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
+    //  <<" range= "<<fRange <<" cMinSt="<<currentMinStep
+    //  << " limit= " << x <<G4endl;
+  }
+  //G4cout<<GetProcessName()<<": e= "<<preStepKinEnergy
   //  <<" stepLimit= "<<x<<G4endl;
   return x;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4VEnergyLossProcess::GetMeanFreePath(
-                             const G4Track& track,
-                             G4double,
-                             G4ForceCondition* condition)
-
-{
-  *condition = NotForced;
-  return MeanFreePath(track);
 }
 
@@ -673 +820 @@
   *condition = NotForced;
   G4double x = DBL_MAX;
+
+  // initialisation of material, mass, charge, model at the beginning of the step
+  DefineMaterial(track.GetMaterialCutsCouple());
+
+  const G4ParticleDefinition* currPart = track.GetDefinition();
+  if(theGenericIon == particle) {
+    massRatio = proton_mass_c2/currPart->GetPDGMass();
+  }  
+  preStepKinEnergy    = track.GetKineticEnergy();
+  preStepScaledEnergy = preStepKinEnergy*massRatio;
+  SelectModel(preStepScaledEnergy);
+
+  if(isIon) {
+    chargeSqRatio = 
+      currentModel->GetChargeSquareRatio(currPart,currentMaterial,preStepKinEnergy);
+    reduceFactor  = 1.0/(chargeSqRatio*massRatio);
+  }
+  //G4cout << "q2= " << chargeSqRatio << " massRatio= " << massRatio << G4endl; 
+  // initialisation for sampling of the interaction length 
   if(previousStepSize <= DBL_MIN) theNumberOfInteractionLengthLeft = -1.0;
-  InitialiseStep(track);
-
+  if(theNumberOfInteractionLengthLeft < 0.0) mfpKinEnergy = DBL_MAX;
+
+  // compute mean free path
   if(preStepScaledEnergy < mfpKinEnergy) {
     if (integral) ComputeLambdaForScaledEnergy(preStepScaledEnergy);
@@ -686 +853 @@
     if (theNumberOfInteractionLengthLeft < 0.0) {
       // beggining of tracking (or just after DoIt of this process)
+      //G4cout<<"G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength Reset"<<G4endl;
       ResetNumberOfInteractionLengthLeft();
     } else if(currentInteractionLength < DBL_MAX) {
@@ -702 +870 @@
       G4cout << "G4VEnergyLossProcess::PostStepGetPhysicalInteractionLength ";
       G4cout << "[ " << GetProcessName() << "]" << G4endl; 
-      G4cout << " for " << particle->GetParticleName() 
+      G4cout << " for " << currPart->GetParticleName() 
              << " in Material  " <<  currentMaterial->GetName()
              << " Ekin(MeV)= " << preStepKinEnergy/MeV 
@@ -710 +878 @@
     }
 #endif
-
     // zero cross section case
   } else {
@@ -738 +905 @@
   if(length <= DBL_MIN) return &fParticleChange;
   G4double eloss  = 0.0;
-
-  /*
+  G4double esecdep = 0.0;
+ 
+  /*  
   if(-1 < verboseLevel) {
     const G4ParticleDefinition* d = track.GetDefinition();
@@ -755 +923 @@
   */
 
+  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
+
   // stopping
   if (length >= fRange) {
+    eloss = preStepKinEnergy;
+    if (useDeexcitation) {
+      if(idxDERegions[currentMaterialIndex]) {
+        currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
+        if(eloss < 0.0) eloss = 0.0;
+      }
+    }
     fParticleChange.SetProposedKineticEnergy(0.0);
-    fParticleChange.ProposeLocalEnergyDeposit(preStepKinEnergy);
+    fParticleChange.ProposeLocalEnergyDeposit(eloss);
     return &fParticleChange;
   }
@@ -772 +949 @@
       GetScaledRangeForScaledEnergy(preStepScaledEnergy) - length/reduceFactor;
     eloss = preStepKinEnergy - ScaledKinEnergyForLoss(x)/massRatio;
+   
     /*
     if(-1 < verboseLevel) 
@@ -781 +959 @@
              << " eloss0(MeV)= "
              << GetDEDXForScaledEnergy(preStepScaledEnergy)*length/MeV
+             << " lim(MeV)= " << preStepKinEnergy*linLossLimit/MeV
              << G4endl;
     */
   }
 
-  const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
-  G4VEmModel* currentModel = SelectModel(preStepScaledEnergy);
-  /*    
+  /*   
   G4double eloss0 = eloss;
   if(-1 < verboseLevel ) {
@@ -801 +978 @@
   G4double cut  = (*theCuts)[currentMaterialIndex];
   G4double esec = 0.0;
-  G4double esecdep = 0.0;
 
   // SubCutOff 
@@ -807 +983 @@
     if(idxSCoffRegions[currentMaterialIndex]) {
 
-      G4double currentMinSafety = 0.0; 
+      G4bool yes = false;
       G4StepPoint* prePoint  = step.GetPreStepPoint();
-      G4StepPoint* postPoint = step.GetPostStepPoint();
-      G4double preSafety  = prePoint->GetSafety();
-      G4double postSafety = preSafety - length; 
-      G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
-
-      // recompute safety
-      if(prePoint->GetStepStatus() != fGeomBoundary &&
-         postPoint->GetStepStatus() != fGeomBoundary) {
-
-        /*
-        //      G4bool yes = (track.GetTrackID() == 5512);
-        G4bool yes = false;
-        if(yes)
-          G4cout << "G4VEnergyLoss: presafety= " << preSafety
-                 << " rcut= " << rcut << "  length= " << length 
-                 << " dir " << track.GetMomentumDirection()
-                 << G4endl;
-        */
-
-        if(preSafety < rcut) 
+
+      // Check boundary
+      if(prePoint->GetStepStatus() == fGeomBoundary) yes = true;
+
+      // Check PrePoint
+      else {
+        G4double preSafety  = prePoint->GetSafety();
+        G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
+
+        // recompute presafety
+        if(preSafety < rcut) {
           preSafety = safetyHelper->ComputeSafety(prePoint->GetPosition());
-
-        //if(yes) {
-        //  G4cout << "G4VEnergyLoss: newsafety= " << preSafety << G4endl;
-          //       if(preSafety==0.0 && track.GetTrackID() == 5512 ) exit(1);
-        //}
-        if(postSafety < rcut) 
-          postSafety = safetyHelper->ComputeSafety(postPoint->GetPosition());
-        /*      
-          if(-1 < verboseLevel) 
-          G4cout << "Subcutoff: presafety(mm)= " << preSafety/mm
-                 << " postsafety(mm)= " << postSafety/mm
-                 << " rcut(mm)= " << rcut/mm 
-                 << G4endl;
-        */
-        currentMinSafety = std::min(preSafety,postSafety); 
-      }
-
+        }
+
+        if(preSafety < rcut) yes = true;
+
+        // Check PostPoint
+        else {
+          G4double postSafety = preSafety - length; 
+          if(postSafety < rcut) {
+            postSafety = 
+              safetyHelper->ComputeSafety(step.GetPostStepPoint()->GetPosition());
+            if(postSafety < rcut) yes = true;
+          }
+        }
+      }
+  
       // Decide to start subcut sampling
-      if(currentMinSafety < rcut) {
+      if(yes) {
 
         cut = (*theSubCuts)[currentMaterialIndex];
@@ -856 +1021 @@
                                 currentModel,currentMaterialIndex, 
                                 esecdep);
+        // add bremsstrahlung sampling
         /*
         if(nProcesses > 0) {
@@ -876 +1042 @@
             esec += e;
             pParticleChange->AddSecondary(t);
-            //  mom -= t->GetMomentum();
           }      
-          //        fParticleChange.SetProposedMomentum(mom);            
         }
       }
@@ -885 +1049 @@
 
   // Corrections, which cannot be tabulated
-  CorrectionsAlongStep(currentCouple, dynParticle, eloss, length);
+  currentModel->CorrectionsAlongStep(currentCouple, dynParticle, 
+                                     eloss, esecdep, length);
 
   // Sample fluctuations
@@ -895 +1060 @@
       G4double tmax = 
         std::min(currentModel->MaxSecondaryKinEnergy(dynParticle),cut);
+      G4double emean = eloss;
       eloss = fluc->SampleFluctuations(currentMaterial,dynParticle,
-                                       tmax,length,eloss);
-      /*           
+                                       tmax,length,emean);
+      /*                            
       if(-1 < verboseLevel) 
       G4cout << "After fluct: eloss(MeV)= " << eloss/MeV
              << " fluc= " << (eloss-eloss0)/MeV
-             << " currentChargeSquare= " << chargeSquare
+             << " ChargeSqRatio= " << chargeSqRatio
              << " massRatio= " << massRatio
              << " tmax= " << tmax
@@ -911 +1077 @@
   eloss += esecdep;
   if(eloss < 0.0) eloss = 0.0;
+
+  // deexcitation
+  else if (useDeexcitation) {
+    if(idxDERegions[currentMaterialIndex]) {
+      currentModel->SampleDeexcitationAlongStep(currentMaterial, track, eloss);
+      if(eloss < 0.0) eloss = 0.0;
+    }
+  }
 
   // Energy balanse
@@ -917 +1091 @@
     eloss  = preStepKinEnergy - esec;
     finalT = 0.0;
+  } else if(isIon) {
+    fParticleChange.SetProposedCharge(
+      currentModel->GetParticleCharge(track.GetDefinition(),currentMaterial,finalT));
   }
 
@@ -922 +1099 @@
   fParticleChange.ProposeLocalEnergyDeposit(eloss);
 
-  /*
+  /*  
   if(-1 < verboseLevel) {
     G4cout << "Final value eloss(MeV)= " << eloss/MeV
@@ -931 +1108 @@
            << G4endl;
   }
-  */
+  */  
 
   return &fParticleChange;
@@ -943 +1120 @@
        G4VEmModel* model,
        G4int idx,
-       G4double& extraEdep) 
+       G4double& /*extraEdep*/) 
 {
   // Fast check weather subcutoff can work
@@ -952 +1129 @@
   const G4Track* track = step.GetTrack();
   const G4DynamicParticle* dp = track->GetDynamicParticle();
+  G4double e = dp->GetKineticEnergy()*massRatio;
   G4bool b;
-  G4double cross = 
-    chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(dp->GetKineticEnergy(),b));
+  G4double cross = chargeSqRatio*(((*theSubLambdaTable)[idx])->GetValue(e,b));
   G4double length = step.GetStepLength();
 
   // negligible probability to get any interaction
   if(length*cross < perMillion) return;
-  /*    
+  /*      
   if(-1 < verboseLevel) 
     G4cout << "<<< Subcutoff for " << GetProcessName()
@@ -970 +1147 @@
   // Sample subcutoff secondaries
   G4StepPoint* preStepPoint = step.GetPreStepPoint();
+  G4StepPoint* postStepPoint = step.GetPostStepPoint();
   G4ThreeVector prepoint = preStepPoint->GetPosition();
-  G4ThreeVector dr = step.GetPostStepPoint()->GetPosition() - prepoint;
+  G4ThreeVector dr = postStepPoint->GetPosition() - prepoint;
   G4double pretime = preStepPoint->GetGlobalTime();
-  //  G4double dt = length/preStepPoint->GetVelocity();
+  G4double dt = postStepPoint->GetGlobalTime() - pretime;
+  //G4double dt = length/preStepPoint->GetVelocity();
   G4double fragment = 0.0;
 
@@ -992 +1171 @@
 
       G4bool addSec = true;
+      /*
       // do not track very low-energy delta-electrons
       if(theSecondaryRangeTable && (*it)->GetDefinition() == theElectron) {
@@ -1004 +1184 @@
         }
       }
+      */
       if(addSec) {
-        //      G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
-        G4Track* t = new G4Track((*it), pretime, r);
+        G4Track* t = new G4Track((*it), pretime + fragment*dt, r);
+        //G4Track* t = new G4Track((*it), pretime, r);
         t->SetTouchableHandle(track->GetTouchableHandle());
         tracks.push_back(t);
 
-        /*
-          if(-1 < verboseLevel) 
-          G4cout << "New track " << p->GetDefinition()->GetParticleName()
-          << " e(keV)= " << p->GetKineticEnergy()/keV
-          << " fragment= " << fragment
-          << G4endl;
+        /*      
+        if(-1 < verboseLevel) 
+          G4cout << "New track " << t->GetDefinition()->GetParticleName()
+                 << " e(keV)= " << t->GetKineticEnergy()/keV
+                 << " fragment= " << fragment
+                 << G4endl;
         */
       }
@@ -1059 +1240 @@
   }
 
-  G4VEmModel* currentModel = SelectModel(postStepScaledEnergy);
+  SelectModel(postStepScaledEnergy);
+  if(useDeexcitation) {
+    currentModel->SetDeexcitationFlag(idxDERegions[currentMaterialIndex]);
+  }
+
   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
   G4double tcut = (*theCuts)[currentMaterialIndex];
@@ -1094 +1279 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-void G4VEnergyLossProcess::PrintInfoDefinition()
-{
-  if(0 < verboseLevel) {
-    G4cout << G4endl << GetProcessName() << ":   tables are built for  "
-           << particle->GetParticleName()
-           << G4endl
-           << "      dE/dx and range tables from "
-           << G4BestUnit(minKinEnergy,"Energy")
-           << " to " << G4BestUnit(maxKinEnergy,"Energy")
-           << " in " << nBins << " bins." << G4endl
-           << "      Lambda tables from threshold to "
-           << G4BestUnit(maxKinEnergy,"Energy")
-           << " in " << nBins << " bins."
+G4bool G4VEnergyLossProcess::StorePhysicsTable(
+       const G4ParticleDefinition* part, const G4String& directory, 
+       G4bool ascii)
+{
+  G4bool res = true;
+  if ( baseParticle || part != particle ) return res;
+
+  if(!StoreTable(part,theDEDXTable,ascii,directory,"DEDX")) 
+    {res = false;}
+
+  if(!StoreTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr")) 
+    {res = false;}
+
+  if(!StoreTable(part,theDEDXSubTable,ascii,directory,"SubDEDX")) 
+    {res = false;}
+
+  if(!StoreTable(part,theIonisationTable,ascii,directory,"Ionisation")) 
+    {res = false;}
+
+  if(!StoreTable(part,theIonisationSubTable,ascii,directory,"SubIonisation")) 
+    {res = false;}
+
+  if(isIonisation &&
+     !StoreTable(part,theCSDARangeTable,ascii,directory,"CSDARange")) 
+    {res = false;}
+
+  if(isIonisation &&
+     !StoreTable(part,theRangeTableForLoss,ascii,directory,"Range")) 
+    {res = false;}
+  
+  if(isIonisation &&
+     !StoreTable(part,theInverseRangeTable,ascii,directory,"InverseRange")) 
+    {res = false;}
+  
+  if(!StoreTable(part,theLambdaTable,ascii,directory,"Lambda")) 
+    {res = false;}
+
+  if(!StoreTable(part,theSubLambdaTable,ascii,directory,"SubLambda")) 
+    {res = false;}
+
+  if ( res ) {
+    if(0 < verboseLevel) {
+      G4cout << "Physics tables are stored for " << particle->GetParticleName()
+             << " and process " << GetProcessName()
+             << " in the directory <" << directory
+             << "> " << G4endl;
+    }
+  } else {
+    G4cout << "Fail to store Physics Tables for " 
+           << particle->GetParticleName()
+           << " and process " << GetProcessName()
+           << " in the directory <" << directory
+           << "> " << G4endl;
+  }
+  return res;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
+       const G4ParticleDefinition* part, const G4String& directory,
+       G4bool ascii)
+{
+  G4bool res = true;
+  const G4String particleName = part->GetParticleName();
+
+  if(1 < verboseLevel) {
+    G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
+           << particleName << " and process " << GetProcessName()
+           << "; tables_are_built= " << tablesAreBuilt
            << G4endl;
-    PrintInfo();
-    if(theRangeTableForLoss && isIonisation) 
-      G4cout << "      Step function: finalRange(mm)= " << finalRange/mm
-             << ", dRoverRange= " << dRoverRange
-             << ", integral: " << integral
-             << ", fluct: " << lossFluctuationFlag
-             << G4endl;
-    
-    if(theCSDARangeTable && isIonisation) 
-      G4cout << "      CSDA range table up"
-             << " to " << G4BestUnit(maxKinEnergyCSDA,"Energy")
-             << " in " << nBinsCSDA << " bins." << G4endl;
-    
-    if(nSCoffRegions>0) 
-      G4cout << "      Subcutoff sampling in " << nSCoffRegions 
-             << " regions" << G4endl;
-
-    if(2 < verboseLevel) {
-      G4cout << "DEDXTable address= " << theDEDXTable << G4endl;
-      if(theDEDXTable && isIonisation) G4cout << (*theDEDXTable) << G4endl;
-      G4cout << "non restricted DEDXTable address= " 
-             << theDEDXunRestrictedTable << G4endl;
-      if(theDEDXunRestrictedTable && isIonisation) 
-           G4cout << (*theDEDXunRestrictedTable) << G4endl;
-      if(theDEDXSubTable && isIonisation) G4cout << (*theDEDXSubTable) 
-                                                 << G4endl;
-      G4cout << "CSDARangeTable address= " << theCSDARangeTable 
+  }
+  if(particle == part) {
+
+    //    G4bool yes = true;
+    if ( !baseParticle ) {
+
+      G4bool fpi = true;
+      if(!RetrieveTable(part,theDEDXTable,ascii,directory,"DEDX",fpi)) 
+        {fpi = false;}
+
+      if(!RetrieveTable(part,theIonisationTable,ascii,directory,"Ionisation",false)) 
+        {fpi = false;}
+
+      if(!RetrieveTable(part,theRangeTableForLoss,ascii,directory,"Range",fpi)) 
+        {res = false;}
+
+      if(!RetrieveTable(part,theDEDXunRestrictedTable,ascii,directory,"DEDXnr",false)) 
+        {res = false;}
+
+      if(!RetrieveTable(part,theCSDARangeTable,ascii,directory,"CSDARange",false)) 
+        {res = false;}
+
+      if(!RetrieveTable(part,theInverseRangeTable,ascii,directory,"InverseRange",fpi)) 
+        {res = false;}
+
+      if(!RetrieveTable(part,theLambdaTable,ascii,directory,"Lambda",true)) 
+        {res = false;}
+
+      G4bool yes = false;
+      if(nSCoffRegions > 0) {yes = true;}
+
+      if(!RetrieveTable(part,theDEDXSubTable,ascii,directory,"SubDEDX",yes)) 
+        {res = false;}
+
+      if(!RetrieveTable(part,theSubLambdaTable,ascii,directory,"SubLambda",yes)) 
+        {res = false;}
+
+      if(!fpi) yes = false;
+      if(!RetrieveTable(part,theIonisationSubTable,ascii,directory,"SubIonisation",yes)) 
+        {res = false;}
+    }
+  }
+
+  return res;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+G4bool G4VEnergyLossProcess::StoreTable(const G4ParticleDefinition* part, 
+                                        G4PhysicsTable* aTable, G4bool ascii,
+                                        const G4String& directory,
+                                        const G4String& tname)
+{
+  G4bool res = true;
+  if ( aTable ) {
+    const G4String name = GetPhysicsTableFileName(part,directory,tname,ascii);
+    if( !aTable->StorePhysicsTable(name,ascii)) res = false;
+  }
+  return res;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
+
+G4bool G4VEnergyLossProcess::RetrieveTable(const G4ParticleDefinition* part, 
+                                           G4PhysicsTable* aTable, G4bool ascii,
+                                           const G4String& directory,
+                                           const G4String& tname,
+                                           G4bool mandatory)
+{
+  G4bool res = true;
+  G4String filename = GetPhysicsTableFileName(part,directory,tname,ascii);
+  G4bool yes = aTable->ExistPhysicsTable(filename);
+  if(yes) {
+    yes = G4PhysicsTableHelper::RetrievePhysicsTable(aTable,filename,ascii);
+    if((G4LossTableManager::Instance())->SplineFlag()) {
+      size_t n = aTable->length();
+      for(size_t i=0; i<n; i++) {(*aTable)[i]->SetSpline(true);}
+    }
+  }
+  if(yes) {
+    if (0 < verboseLevel) {
+      G4cout << tname << " table for " << part->GetParticleName() 
+             << " is Retrieved from <" << filename << ">"
              << G4endl;
-      if(theCSDARangeTable && isIonisation) G4cout << (*theCSDARangeTable) 
-            << G4endl;
-      G4cout << "RangeTableForLoss address= " << theRangeTableForLoss 
+    }
+  } else {
+    if(mandatory) res = false;
+    if(mandatory || 1 < verboseLevel) {
+      G4cout << tname << " table for " << part->GetParticleName() 
+             << " from file <"
+             << filename << "> is not Retrieved"
              << G4endl;
-      if(theRangeTableForLoss && isIonisation) 
-             G4cout << (*theRangeTableForLoss) << G4endl;
-      G4cout << "InverseRangeTable address= " << theInverseRangeTable 
-             << G4endl;
-      if(theInverseRangeTable && isIonisation) 
-             G4cout << (*theInverseRangeTable) << G4endl;
-      G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
-      if(theLambdaTable && isIonisation) G4cout << (*theLambdaTable) << G4endl;
-      G4cout << "SubLambdaTable address= " << theSubLambdaTable << G4endl;
-      if(theSubLambdaTable && isIonisation) G4cout << (*theSubLambdaTable) 
-             << G4endl;
+    }
+  }
+  return res;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEnergyLossProcess::GetDEDXDispersion(
+                                  const G4MaterialCutsCouple *couple,
+                                  const G4DynamicParticle* dp,
+                                        G4double length)
+{
+  DefineMaterial(couple);
+  G4double ekin = dp->GetKineticEnergy();
+  SelectModel(ekin*massRatio);
+  G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
+  tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
+  G4double d = 0.0;
+  G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
+  if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
+  return d;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEnergyLossProcess::CrossSectionPerVolume(
+         G4double kineticEnergy, const G4MaterialCutsCouple* couple)
+{
+  // Cross section per volume is calculated
+  DefineMaterial(couple);
+  G4double cross = 0.0;
+  G4bool b;
+  if(theLambdaTable) {
+    cross = 
+      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b);
+  } else {
+    SelectModel(kineticEnergy);
+    cross = 
+      currentModel->CrossSectionPerVolume(currentMaterial,
+                                          particle, kineticEnergy,
+                                          (*theCuts)[currentMaterialIndex]);
+  }
+
+  return cross;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEnergyLossProcess::MeanFreePath(const G4Track& track)
+{
+  DefineMaterial(track.GetMaterialCutsCouple());
+  preStepLambda = GetLambdaForScaledEnergy(track.GetKineticEnergy()*massRatio);
+  G4double x = DBL_MAX;
+  if(DBL_MIN < preStepLambda) x = 1.0/preStepLambda;
+  return x;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEnergyLossProcess::ContinuousStepLimit(const G4Track& track, 
+                                                   G4double x, G4double y, 
+                                                   G4double& z)
+{
+  G4GPILSelection sel;
+  return AlongStepGetPhysicalInteractionLength(track, x, y, z, &sel);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEnergyLossProcess::GetMeanFreePath(
+                             const G4Track& track,
+                             G4double,
+                             G4ForceCondition* condition)
+
+{
+  *condition = NotForced;
+  return MeanFreePath(track);
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4double G4VEnergyLossProcess::GetContinuousStepLimit(
+                const G4Track&,
+                G4double, G4double, G4double&)
+{
+  return DBL_MAX;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+
+G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
+                 const G4MaterialCutsCouple* couple, G4double cut)
+{
+  //  G4double cut  = (*theCuts)[couple->GetIndex()];
+  //  G4int nbins = nLambdaBins;
+  G4double tmin = 
+    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
+             minKinEnergy);
+  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
+  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
+  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
+
+  return v;
+}
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
+  
+void G4VEnergyLossProcess::AddCollaborativeProcess(
+            G4VEnergyLossProcess* p)
+{
+  G4bool add = true;
+  if(p->GetProcessName() != "eBrem") add = false;
+  if(add && nProcesses > 0) {
+    for(G4int i=0; i<nProcesses; i++) {
+      if(p == scProcesses[i]) {
+        add = false;
+        break;
+      }
+    }
+  }
+  if(add) {
+    scProcesses.push_back(p);
+    nProcesses++;
+    if (1 < verboseLevel) { 
+      G4cout << "### The process " << p->GetProcessName() 
+             << " is added to the list of collaborative processes of "
+             << GetProcessName() << G4endl; 
     }
   }
@@ -1182 +1600 @@
   } else if(fSubRestricted == tType) {    
     theDEDXSubTable = p;
-  } else if(fIonisation == tType && theIonisationTable != p) {    
+  } else if(fIsIonisation == tType && theIonisationTable != p) {    
     if(theIonisationTable) theIonisationTable->clearAndDestroy();
     theIonisationTable = p;
-  } else if(fSubIonisation == tType && theIonisationSubTable != p) {    
+  } else if(fIsSubIonisation == tType && theIonisationSubTable != p) {    
     if(theIonisationSubTable) theIonisationSubTable->clearAndDestroy();
     theIonisationSubTable = p;
@@ -1319 +1737 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4PhysicsVector* G4VEnergyLossProcess::LambdaPhysicsVector(
-                 const G4MaterialCutsCouple* couple, G4double cut)
-{
-  //  G4double cut  = (*theCuts)[couple->GetIndex()];
-  //  G4int nbins = nLambdaBins;
-  G4double tmin = 
-    std::max(MinPrimaryEnergy(particle, couple->GetMaterial(), cut),
-             minKinEnergy);
-  if(tmin >= maxKinEnergy) tmin = 0.5*maxKinEnergy;
-  G4PhysicsVector* v = new G4PhysicsLogVector(tmin, maxKinEnergy, nBins);
-  return v;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4VEnergyLossProcess::MicroscopicCrossSection(
-         G4double kineticEnergy, const G4MaterialCutsCouple* couple)
-{
-  // Cross section per atom is calculated
-  DefineMaterial(couple);
-  G4double cross = 0.0;
-  G4bool b;
-  if(theLambdaTable) 
-    cross = 
-      ((*theLambdaTable)[currentMaterialIndex])->GetValue(kineticEnergy, b)/
-      currentMaterial->GetTotNbOfAtomsPerVolume();
-
-  return cross;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4bool G4VEnergyLossProcess::StorePhysicsTable(
-       const G4ParticleDefinition* part, const G4String& directory, 
-       G4bool ascii)
-{
-  G4bool res = true;
-  if ( baseParticle || part != particle ) return res;
-
-  if ( theDEDXTable ) {
-    const G4String name = GetPhysicsTableFileName(part,directory,"DEDX",ascii);
-    if( !theDEDXTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theDEDXunRestrictedTable ) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"DEDXnr",ascii);
-    if( !theDEDXTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theDEDXSubTable ) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"SubDEDX",ascii);
-    if( !theDEDXSubTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theIonisationTable ) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"Ionisation",ascii);
-    if( !theIonisationTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theIonisationSubTable ) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"SubIonisation",ascii);
-    if( !theIonisationSubTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theCSDARangeTable && isIonisation ) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"CSDARange",ascii);
-    if( !theCSDARangeTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theRangeTableForLoss && isIonisation ) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"Range",ascii);
-    if( !theRangeTableForLoss->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theInverseRangeTable && isIonisation ) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"InverseRange",ascii);
-    if( !theInverseRangeTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theLambdaTable  && isIonisation) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"Lambda",ascii);
-    if( !theLambdaTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-
-  if ( theSubLambdaTable  && isIonisation) {
-    const G4String name = 
-      GetPhysicsTableFileName(part,directory,"SubLambda",ascii);
-    if( !theSubLambdaTable->StorePhysicsTable(name,ascii)) res = false;
-  }
-  if ( res ) {
-    if(0 < verboseLevel) {
-      G4cout << "Physics tables are stored for " << particle->GetParticleName()
-             << " and process " << GetProcessName()
-             << " in the directory <" << directory
-             << "> " << G4endl;
-    }
-  } else {
-    G4cout << "Fail to store Physics Tables for " 
-           << particle->GetParticleName()
-           << " and process " << GetProcessName()
-           << " in the directory <" << directory
-           << "> " << G4endl;
-  }
-  return res;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.....
-
-G4bool G4VEnergyLossProcess::RetrievePhysicsTable(
-       const G4ParticleDefinition* part, const G4String& directory,
-       G4bool ascii)
-{
-  G4bool res = true;
-  const G4String particleName = part->GetParticleName();
-
-  if(1 < verboseLevel) {
-    G4cout << "G4VEnergyLossProcess::RetrievePhysicsTable() for "
-           << particleName << " and process " << GetProcessName()
-           << "; tables_are_built= " << tablesAreBuilt
-           << G4endl;
-  }
-  if(particle == part) {
-
-    G4bool yes = true;
-    G4bool fpi = true;
-    if ( !baseParticle ) {
-      G4String filename;
-
-      filename = GetPhysicsTableFileName(part,directory,"DEDX",ascii);
-      yes = theDEDXTable->ExistPhysicsTable(filename);
-      if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theDEDXTable,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "DEDX table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        fpi = false;
-        if (1 < verboseLevel) {
-          G4cout << "DEDX table for " << particleName << " from file <"
-                 << filename << "> is not Retrieved"
-                 << G4endl;
-        }
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"Range",ascii);
-      yes = theRangeTableForLoss->ExistPhysicsTable(filename);
-      if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theRangeTableForLoss,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "Range table for loss for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        if(fpi) {
-          res = false;
-          G4cout << "Range table for loss for " << particleName 
-                 << " from file <"
-                 << filename << "> is not Retrieved"
-                 << G4endl;
-        }
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"DEDXnr",ascii);
-      yes = theDEDXunRestrictedTable->ExistPhysicsTable(filename);
-      if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theDEDXunRestrictedTable,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "Non-restricted DEDX table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        if (1 < verboseLevel) {
-          G4cout << "Non-restricted DEDX table for " << particleName 
-                 << " from file <"
-                 << filename << "> is not Retrieved"
-                 << G4endl;
-        }
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"CSDARange",ascii);
-      yes = theCSDARangeTable->ExistPhysicsTable(filename);
-      if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theCSDARangeTable,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "CSDA Range table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        G4cout << "CSDA Range table for loss for " << particleName 
-               << " does not exist"
-               << G4endl;
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"InverseRange",ascii);
-      yes = theInverseRangeTable->ExistPhysicsTable(filename);
-      if(yes)  yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                     theInverseRangeTable,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "InverseRange table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        if(fpi) {
-          res = false;
-          G4cout << "InverseRange table for " << particleName 
-                 << " from file <"
-                 << filename << "> is not Retrieved"
-                 << G4endl;
-
-        }
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
-      yes = theLambdaTable->ExistPhysicsTable(filename);
-      if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theLambdaTable,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "Lambda table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        if(fpi) {
-          res = false;
-          G4cout << "Lambda table for " << particleName << " from file <"
-                 << filename << "> is not Retrieved"
-                 << G4endl;
-        }
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"SubDEDX",ascii);
-      yes = theDEDXSubTable->ExistPhysicsTable(filename);
-      if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theDEDXSubTable,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "SubDEDX table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        if(nSCoffRegions) {
-          res=false;
-          G4cout << "SubDEDX table for " << particleName << " from file <"
-                 << filename << "> is not Retrieved"
-                 << G4endl;
-        }
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"SubLambda",ascii);
-      yes = theSubLambdaTable->ExistPhysicsTable(filename);
-      if(yes) yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theSubLambdaTable,filename,ascii);
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "SubLambda table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } else {
-        if(nSCoffRegions) {
-          res=false;
-          G4cout << "SubLambda table for " << particleName << " from file <"
-                 << filename << "> is not Retrieved"
-                 << G4endl;
-        }
-      }
-
-      filename = GetPhysicsTableFileName(part,directory,"Ionisation",ascii);
-      yes = theIonisationTable->ExistPhysicsTable(filename);
-      if(yes) {
-        theIonisationTable = 
-          G4PhysicsTableHelper::PreparePhysicsTable(theIonisationTable);
-        
-        yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theIonisationTable,filename,ascii);
-      }
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "Ionisation table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } 
-
-      filename = GetPhysicsTableFileName(part,directory,"SubIonisation",ascii);
-      yes = theIonisationSubTable->ExistPhysicsTable(filename);
-      if(yes) {
-        theIonisationSubTable = 
-          G4PhysicsTableHelper::PreparePhysicsTable(theIonisationSubTable);
-        yes = G4PhysicsTableHelper::RetrievePhysicsTable(
-                    theIonisationSubTable,filename,ascii);
-      }
-      if(yes) {
-        if (0 < verboseLevel) {
-          G4cout << "SubIonisation table for " << particleName 
-                 << " is Retrieved from <"
-                 << filename << ">"
-                 << G4endl;
-        }
-      } 
-    }
-  }
-
-  return res;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-  
-void G4VEnergyLossProcess::AddCollaborativeProcess(
-            G4VEnergyLossProcess* p)
-{
-  G4bool add = true;
-  if(nProcesses > 0) {
-    for(G4int i=0; i<nProcesses; i++) {
-      if(p == scProcesses[i]) {
-        add = false;
-        break;
-      }
-    }
-  }
-  if(add) {
-    scProcesses.push_back(p);
-    nProcesses++;
-    if (0 < verboseLevel) 
-      G4cout << "### The process " << p->GetProcessName() 
-             << " is added to the list of collaborative processes of "
-             << GetProcessName() << G4endl; 
-  }
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-G4double G4VEnergyLossProcess::GetDEDXDispersion(
-                                  const G4MaterialCutsCouple *couple,
-                                  const G4DynamicParticle* dp,
-                                        G4double length)
-{
-  DefineMaterial(couple);
-  G4double ekin = dp->GetKineticEnergy();
-  G4VEmModel* currentModel = SelectModel(ekin*massRatio);
-  G4double tmax = currentModel->MaxSecondaryKinEnergy(dp);
-  tmax = std::min(tmax,(*theCuts)[currentMaterialIndex]);
-  G4double d = 0.0;
-  G4VEmFluctuationModel* fm = currentModel->GetModelOfFluctuations();
-  if(fm) d = fm->Dispersion(currentMaterial,dp,tmax,length);
-  return d;
-}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-
-void G4VEnergyLossProcess::ActivateDeexcitation(G4bool, const G4Region*)
-{}
-
-//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
-

trunk/source/processes/electromagnetic/utils/src/G4VMultipleScattering.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4VMultipleScattering.cc,v 1.47 2007/11/09 11:35:54 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4VMultipleScattering.cc,v 1.60 2008/11/20 20:32:40 vnivanch Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -55 +55 @@
 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
+// 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko)
 //
 // Class Description:
@@ -85 +86 @@
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 
-G4VMultipleScattering::G4VMultipleScattering(const G4String& name, G4ProcessType type):
-                 G4VContinuousDiscreteProcess(name, type),
+G4VMultipleScattering::G4VMultipleScattering(const G4String& name, 
+                                             G4ProcessType type):
+  G4VContinuousDiscreteProcess(name, type),
+  buildLambdaTable(true),
   theLambdaTable(0),
   firstParticle(0),
-  currentParticle(0),
-  currentCouple(0),
-  nBins(120),
   stepLimit(fUseSafety),
-  skin(0.0),
+  skin(3.0),
   facrange(0.02),
   facgeom(2.5),
   latDisplasment(true),
-  buildLambdaTable(true)
-{
+  currentParticle(0),
+  currentCouple(0)
+{
+  SetVerboseLevel(1);
+  SetProcessSubType(fMultipleScattering);
+
+  // Size of tables assuming spline
   minKinEnergy = 0.1*keV;
   maxKinEnergy = 100.0*TeV;
-  SetVerboseLevel(1);
+  nBins        = 84;
+
+  // default limit on polar angle
+  polarAngleLimit = 0.0;
 
   pParticleChange = &fParticleChange;
@@ -114 +122 @@
 G4VMultipleScattering::~G4VMultipleScattering()
 {
-  if(1 < verboseLevel) 
+  if(1 < verboseLevel) {
     G4cout << "G4VMultipleScattering destruct " << GetProcessName() 
            << G4endl;
+  }
   delete modelManager;
   if (theLambdaTable) {
@@ -131 +140 @@
   G4String num = part.GetParticleName();
   if(1 < verboseLevel) {
-    //    G4cout << "========================================================" << G4endl;
     G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
            << GetProcessName()
@@ -148 +156 @@
       if (theLambdaTable->GetFlag(i)) {
         // create physics vector and fill it
-        const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
+        const G4MaterialCutsCouple* couple = 
+          theCoupleTable->GetMaterialCutsCouple(i);
         G4PhysicsVector* aVector = PhysicsVector(couple);
         modelManager->FillLambdaVector(aVector, couple, false);
@@ -162 +171 @@
   }
   if(verboseLevel>0 && ( num == "e-" || num == "mu+" ||  
-                         num == "proton" || num == "pi-" || num == "GenericIon")) {
+                         num == "proton" || num == "pi-" || 
+                         num == "GenericIon")) {
     PrintInfoDefinition();
     if(2 < verboseLevel && theLambdaTable) G4cout << *theLambdaTable << G4endl;
@@ -182 +192 @@
     currentCouple = 0;
     if(part.GetParticleType() == "nucleus" && 
-       part.GetParticleSubType() == "generic") 
-         firstParticle = G4GenericIon::GenericIon();
-    else firstParticle = &part; 
+       part.GetParticleSubType() == "generic") {
+      firstParticle = G4GenericIon::GenericIon();
+    } else {
+      firstParticle = &part;
+    } 
+
     currentParticle = &part;
   }
 
   if(1 < verboseLevel) {
-    //    G4cout << "========================================================" << G4endl;
     G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
            << GetProcessName()
@@ -217 +229 @@
 {
   if (0 < verboseLevel) {
-    G4cout << G4endl << GetProcessName() << ":  Model variant of multiple scattering "
-           << "for " << firstParticle->GetParticleName()
+    G4cout << G4endl << GetProcessName() 
+           << ":   for " << firstParticle->GetParticleName()
+           << "    SubType= " << GetProcessSubType() 
            << G4endl;
     if (theLambdaTable) {
@@ -225 +238 @@
              << " to "
              << G4BestUnit(MaxKinEnergy(),"Energy")
-             << " in " << nBins << " bins."
+             << " in " << nBins << " bins, spline: " 
+             << (G4LossTableManager::Instance())->SplineFlag()
              << G4endl;
     }
-    G4cout << "      LateralDisplacementFlag=  " << latDisplasment
-           << "   Skin= " << skin << G4endl; 
     PrintInfo();
+    modelManager->DumpModelList(verboseLevel);
     if (2 < verboseLevel) {
       G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
@@ -242 +255 @@
 G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength(
                              const G4Track& track,
-                             G4double previousStepSize,
+                             G4double,
                              G4double currentMinimalStep,
                              G4double& currentSafety,
@@ -249 +262 @@
   // get Step limit proposed by the process
   valueGPILSelectionMSC = NotCandidateForSelection;
-  G4double steplength = GetMscContinuousStepLimit(track,previousStepSize,
-                                              currentMinimalStep,currentSafety);
+  G4double steplength = GetMscContinuousStepLimit(track,
+                                                  track.GetKineticEnergy(),
+                                                  currentMinimalStep,
+                                                  currentSafety);
   // G4cout << "StepLimit= " << steplength << G4endl;
   // set return value for G4GPILSelection
@@ -315 +330 @@
   if( couple->IsUsed() ) nbins = nBins;
   G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nbins);
+  v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
   return v;
 }
@@ -372 +388 @@
                << G4endl;
     }
+    if((G4LossTableManager::Instance())->SplineFlag()) {
+      size_t n = theLambdaTable->length();
+      for(size_t i=0; i<n; i++) {(* theLambdaTable)[i]->SetSpline(true);}
+    }
   } else {
     if (1 < verboseLevel) {

trunk/source/processes/electromagnetic/utils/src/G4ionEffectiveCharge.cc

@@ -24 +24 @@
 // ********************************************************************
 //
-// $Id: G4ionEffectiveCharge.cc,v 1.17.2.1 2008/04/22 15:28:13 vnivanch Exp $
-// GEANT4 tag $Name: geant4-09-01-patch-02 $
+// $Id: G4ionEffectiveCharge.cc,v 1.24 2008/12/18 13:01:46 gunter Exp $
+// GEANT4 tag $Name: geant4-09-02-ref-02 $
 //
 // -------------------------------------------------------------------
@@ -54 +54 @@
 #include "G4ionEffectiveCharge.hh"
 #include "G4UnitsTable.hh"
-#include "G4ParticleDefinition.hh"
 #include "G4Material.hh"
+#include "G4NistManager.hh"
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
@@ -62 +62 @@
 {
   chargeCorrection = 1.0;
-  energyHighLimit  = 10.0*MeV;
+  energyHighLimit  = 20.0*MeV;
   energyLowLimit   = 1.0*keV;
   energyBohr       = 25.*keV;
   massFactor       = amu_c2/(proton_mass_c2*keV);
-  minCharge        = 0.1;
+  minCharge        = 1.0;
+  lastPart         = 0;
+  lastMat          = 0;
+  lastKinEnergy    = 0.0;
+  effCharge        = eplus;
+  nist = G4NistManager::Instance();
 }
 
@@ -78 +83 @@
 G4double G4ionEffectiveCharge::EffectiveCharge(const G4ParticleDefinition* p,
                                                const G4Material* material,
-                                                     G4double kineticEnergy)
+                                               G4double kineticEnergy)
 {
+  if(p == lastPart && material == lastMat && kineticEnergy == lastKinEnergy)
+    return effCharge;
+
+  lastPart      = p;
+  lastMat       = material;
+  lastKinEnergy = kineticEnergy;
+
   G4double mass   = p->GetPDGMass();
   G4double charge = p->GetPDGCharge();
@@ -85 +97 @@
 
   chargeCorrection = 1.0;
+  effCharge = charge;
 
   // The aproximation of ion effective charge from:
@@ -94 +107 @@
   if( reducedEnergy > Zi*energyHighLimit || Zi < 1.5 || !material) return charge;
 
-  static G4double c[6] = {0.2865,  0.1266, -0.001429,
-                          0.02402,-0.01135, 0.001475} ;
-
   G4double z    = material->GetIonisation()->GetZeffective();
-  reducedEnergy = std::max(reducedEnergy,energyLowLimit);
-  G4double q;
+  //  reducedEnergy = std::max(reducedEnergy,energyLowLimit);
 
   // Helium ion case
   if( Zi < 2.5 ) {
+
+    static G4double c[6] = {0.2865,  0.1266, -0.001429,
+                            0.02402,-0.01135, 0.001475} ;
 
     G4double Q = std::max(0.0,std::log(reducedEnergy*massFactor));
@@ -120 +132 @@
     if(tq2 < 0.2) tt *= (1.0 - tq2 + 0.5*tq2*tq2);
     else          tt *= std::exp(-tq2);
-    q = (1.0 + tt) * std::sqrt(ex);
+
+    effCharge = charge*(1.0 + tt) * std::sqrt(ex);
 
     // Heavy ion case
   } else {
     
-    G4double z23  = std::pow(z, 0.666666);
-    G4double zi13 = std::pow(Zi, 0.333333);
+    G4double y;
+    //    = nist->GetZ13(z);
+    //G4double z23  = y*y;
+    G4double zi13 = nist->GetZ13(Zi);
     G4double zi23 = zi13*zi13;
-    G4double e = std::max(reducedEnergy,energyBohr/z23);
+    //    G4double e = std::max(reducedEnergy,energyBohr/z23);
+    //G4double e = reducedEnergy;
 
     // v1 is ion velocity in vF unit
     G4double eF   = material->GetIonisation()->GetFermiEnergy();
-    G4double v1sq = e/eF;
+    G4double v1sq = reducedEnergy/eF;
     G4double vFsq = eF/energyBohr;
-    G4double vF   = std::sqrt(vFsq);
-
-    G4double y ;
+    G4double vF   = std::sqrt(eF/energyBohr);
 
     // Faster than Fermi velocity
@@ -147 +161 @@
     }
 
+    G4double q;
     G4double y3 = std::pow(y, 0.3) ;
-    //    G4cout << "y= " << y << " y3= " << y3 << " v1= " << v1 << " vF= " << vF << G4endl; 
+    // G4cout<<"y= "<<y<<" y3= "<<y3<<" v1= "<<v1<<" vF= "<<vF<<G4endl; 
     q = 1.0 - std::exp( 0.803*y3 - 1.3167*y3*y3 - 0.38157*y - 0.008983*y*y ) ;
-
-    //    q = 1.0 - std::exp(-0.95*std::sqrt(reducedEnergy/energyBohr)/zi23);
+    
+    //y *= 0.77;
+    //y *= (0.75 + 0.52/Zi);
+
+    //if( y < 0.2 ) q = y*(1.0 - 0.5*y);
+    //else          q = 1.0 - std::exp(-y);
 
     G4double qmin = minCharge/Zi;
-
     if(q < qmin) q = qmin;
+  
+    effCharge = q*charge;
+
+    /*
+    G4double x1 = 1.0*effCharge*(1.0 - 0.132*std::log(y))/(y*std::sqrt(z));
+    G4double x2 = 0.1*effCharge*effCharge*energyBohr/reducedEnergy;
+
+    chargeCorrection = 1.0 + x1 - x2;
+
+    G4cout << "x1= "<<x1<<" x2= "<< x2<<" corr= "<<chargeCorrection<<G4endl;
+    */
     
     G4double tq = 7.6 - std::log(reducedEnergy/keV);
@@ -170 +199 @@
 
     G4double lambda = 10.0 * vF / (zi13 * (6.0 + q));
-    if(q < 0.2) lambda *= (1.0 - 0.666666*q - q*q/9.0);
+    if(q < 0.2) lambda *= (1.0 - 0.66666667*q - q*q/9.0);
     else        lambda *= std::pow(1.0-q, 0.666666);
 
@@ -180 +209 @@
 
     chargeCorrection = sq * (1.0 + xx);
+    
   }
   //  G4cout << "G4ionEffectiveCharge: charge= " << charge << " q= " << q 
   //         << " chargeCor= " << chargeCorrection 
   //       << " e(MeV)= " << kineticEnergy/MeV << G4endl;
-  return q*charge;
+  return effCharge;
 }