COMP/CMSDIST/geant4.9.5.p01-fix-cputime-regression.patch

diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh geant4.9.5.p01/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh
--- geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh       2012-07-23 08:38:19.000000000 +0200
+++ geant4.9.5.p01/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh    2012-07-23 09:14:23.000000000 +0200
@@ -58,6 +58,8 @@
 
 enum { TSIZE=41, NPARTS=35, NELAB=17, NCNLW=15, NFISS=21 };
 
+class G4Pow;
+
 class G4HadronCrossSections
 {
   public:
@@ -100,6 +102,8 @@ class G4HadronCrossSections
 
     static G4HadronCrossSections* theInstance;
 
+    G4Pow* g4pow;
+
     G4double sigelastic;
     G4double siginelastic;
     G4ParticleDefinition* prevParticleDefinition;
diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc geant4.9.5.p01/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc
--- geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc   2012-07-23 08:38:19.000000000 +0200
+++ geant4.9.5.p01/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc        2012-07-23 09:14:05.000000000 +0200
@@ -55,7 +55,7 @@
 #include "G4HadronCrossSections.hh"
 #include "G4ios.hh"
 #include "G4HadronicException.hh"
-
+#include "G4Pow.hh"
  
 // Initialize static pointer for singleton instance
 G4HadronCrossSections* G4HadronCrossSections::theInstance = 0;
@@ -1217,7 +1217,10 @@ G4bool G4HadronCrossSections::correctIne
 G4HadronCrossSections::G4HadronCrossSections()
   : prevParticleDefinition(0), prevElement(0), prevZZ(0), prevAA(0), 
     prevKineticEnergy(DBL_MAX), lastEkx(0.), lastEkxPower(0.), verboseLevel(0)
-{}
+{
+  siginelastic = sigelastic = 0.0;
+  g4pow = G4Pow::GetInstance();
+}
 
 G4HadronCrossSections::~G4HadronCrossSections() 
 {}
@@ -1243,7 +1246,12 @@ G4HadronCrossSections::GetInelasticCross
   if (particle->GetDefinition() != prevParticleDefinition ||
       particle->GetKineticEnergy() != prevKineticEnergy ||
       ZZ != prevZZ || AA != prevAA)
-    CalcScatteringCrossSections(particle, ZZ, AA);
+    { CalcScatteringCrossSections(particle, ZZ, AA); }
+  if (verboseLevel > 1) {
+    G4cout << "        x-section inelastic(mb)= " 
+          << siginelastic/millibarn << G4endl; 
+  }
+
   return siginelastic;
 }
 
@@ -1255,8 +1263,12 @@ G4HadronCrossSections::GetElasticCrossSe
   if (particle->GetDefinition() != prevParticleDefinition ||
       particle->GetKineticEnergy() != prevKineticEnergy ||
       ZZ != prevZZ || AA != prevAA)
-    CalcScatteringCrossSections(particle, ZZ, AA);
+    { CalcScatteringCrossSections(particle, ZZ, AA); }
 
+  if (verboseLevel > 1) {
+    G4cout << "        x-section elastic(mb)= " 
+          << sigelastic/millibarn << G4endl; 
+  }
   return sigelastic;
 }
 
@@ -1279,33 +1291,37 @@ G4HadronCrossSections::CalcScatteringCro
    G4double xspiel, xspiin;
 
    G4int ipart = GetParticleCode(aParticle);
-   G4double a = AA;
-   G4double z = ZZ;
+   prevAA = AA;
+   prevZZ = ZZ;
+   prevParticleDefinition = aParticle->GetDefinition(); 
+   prevKineticEnergy = aParticle->GetKineticEnergy(); 
 
    if (verboseLevel > 1) {
-      G4cout << "G4HadronCrossSections: a=" << a << G4endl;
-      G4cout << "G4HadronCrossSections: z=" << z << G4endl;
+     G4cout << "G4HadronCrossSections::ScatteringCrossSections: " 
+           << aParticle->GetDefinition()->GetParticleName()
+           << " E(MeV)= " << aParticle->GetKineticEnergy()/MeV 
+           << " Z= " << ZZ << " A= " << AA << G4endl;
    }
 
-// Ions...
+   // Ions...
 
    if (ipart >= 30 && ipart <= 32) {
 
       G4double apart=0;
-      if (ipart == 30) apart = std::pow(2., 1./3.);
-      else if (ipart == 31) apart = std::pow(3., 1./3.);
-      else if (ipart == 32) apart = std::pow(4., 1./3.);
+      if (ipart == 30) apart = g4pow->Z13(2);
+      else if (ipart == 31) apart = g4pow->Z13(3);
+      else if (ipart == 32) apart = g4pow->Z13(4);
 
-      G4double term = apart + std::pow(a, 1./3.);
+      G4double term = apart + g4pow->Z13(AA);
       sigin = 49.*term*term;
-   // Convert cross section from mb to default units
+      // Convert cross section from mb to default units
       siginelastic = sigin*millibarn;
-      if(aParticle->GetKineticEnergy() < 6*MeV) siginelastic = 0;
+      if(prevKineticEnergy < 6*MeV) siginelastic = 0;
       sigelastic = 0.;
       return;
    }
 
-   G4double ek = aParticle->GetKineticEnergy()/GeV;
+   G4double ek = prevKineticEnergy/GeV;
 
 // Low energy neutrons...
 
@@ -1349,7 +1365,7 @@ G4HadronCrossSections::CalcScatteringCro
       G4int ja2 = NCNLW - 1;
       do {
          G4int midBin = (ja1 + ja2)/2;
-         if (a < cnlwat[midBin])
+         if (AA < cnlwat[midBin])
            ja2 = midBin;
          else
            ja1 = midBin;
@@ -1365,7 +1381,7 @@ G4HadronCrossSections::CalcScatteringCro
       dy = cnlwel[ja2][je1] - cnlwel[ja1][je1];
       rca = dy/dnlwat;
       b = cnlwel[ja1][je1] - rce*elab[je1] - rca*cnlwat[ja1];
-      sigelastic = rce*ek + rca*a + b;
+      sigelastic = rce*ek + rca*AA + b;
 // Inelastic cross section:
       // E interpolation or extrapolation at JA1
       dy = cnlwin[ja1][je2] - cnlwin[ja1][je1];
@@ -1374,7 +1390,7 @@ G4HadronCrossSections::CalcScatteringCro
       dy = cnlwin[ja2][je1] - cnlwin[ja1][je1];
       rca = dy/dnlwat;
       b = cnlwin[ja1][je1] - rce*elab[je1] - rca*cnlwat[ja1];
-      siginelastic = rce*ek + rca*a + b;
+      siginelastic = rce*ek + rca*AA + b;
    // Convert cross sections from mb to default units
       sigelastic = sigelastic*millibarn;
       siginelastic = siginelastic*millibarn;
@@ -1446,14 +1462,14 @@ G4HadronCrossSections::CalcScatteringCro
 
 // A-dependence from parameterization...
 
-   if (a >= 1.5) {
+   if (AA >= 2) {
 
       crel = 1.;
       crin = 1.;
 
       G4int i = 3;
-      if (a < 50.) i = 2;
-      if (a > 100.) i = 4;
+      if (AA < 50) i = 2;
+      else if (AA > 100.) i = 4;
       i = i - 1;      // For array indexing
 
 // Protons and neutrons
@@ -1474,7 +1490,7 @@ G4HadronCrossSections::CalcScatteringCro
          if (correctInelasticNearZero && je1 == 0 && i <= 3) {
             G4double m0 = aParticle->GetMass()/GeV;
             G4double T = std::sqrt(m0*m0 + p*p) - m0;
-            G4double dx = std::sqrt(m0*m0 + plab[1]*plab[1]) - m0;
+            dx = std::sqrt(m0*m0 + plab[1]*plab[1]) - m0;
             rc = dy/dx;
             xsecin = rc*T + b;
          }
@@ -1490,7 +1506,7 @@ G4HadronCrossSections::CalcScatteringCro
 // Calculate correction factors (crel, crin) from values
 // on Al, Cu, Pb.  Note that data is only available for pions and protons.
          G4double wgch = 0.5;
-         if (a < 20.) wgch = 0.5 + 0.5*std::exp(-(a - 1.));
+         if (AA < 20) { wgch = 0.5 + 0.5*std::exp(-G4double(AA - 1)); }
          sigel = wgch*sigel + (1. - wgch)*xsecel;
          sigin = wgch*sigin + (1. - wgch)*xsecin;
          
@@ -1517,8 +1533,8 @@ G4HadronCrossSections::CalcScatteringCro
       }
 
 // Apply correction factors
-      sigin = crin*(sigin + sigel)*std::pow(a, alph);
-      sigel = crel*0.36*sigel*std::pow(a, 1.17);
+      sigin = crin*(sigin + sigel)*g4pow->powZ(AA, alph);
+      sigel = crel*0.36*sigel*g4pow->powZ(AA, 1.17);
       sigel = sigel*partel[ipart1];
       sigin = sigin*partin[ipart1];
    }
@@ -1591,7 +1607,7 @@ G4HadronCrossSections::GetFissionCrossSe
 
   G4double z43ba;
   if (j == 4) {
-    z43ba = std::pow(G4double(ZZ), 4./3.)/G4double(AA);
+    z43ba = ZZ*g4pow->Z13(ZZ)/G4double(AA);
     z43ba = std::max(-67. + 38.7*z43ba, 0.);
   } else {
     z43ba = 1.;
diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/management/src/G4HadronicProcess.cc geant4.9.5.p01/source/processes/hadronic/management/src/G4HadronicProcess.cc
--- geant4.9.5.p01.orig/source/processes/hadronic/management/src/G4HadronicProcess.cc   2012-07-23 08:38:33.000000000 +0200
+++ geant4.9.5.p01/source/processes/hadronic/management/src/G4HadronicProcess.cc        2012-07-23 09:15:11.000000000 +0200
@@ -213,11 +213,13 @@ G4HadronicProcess::PostStepDoIt(const G4
                ed);
   }
 
-  if (GetElementCrossSection(aParticle, anElement, aMaterial) <= 0.0) {
-    // No interaction
-    //theTotalResult->Clear();
-    return theTotalResult;
-  }    
+  if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) {
+    if (GetElementCrossSection(aParticle, anElement, aMaterial) <= 0.0) {
+      // No interaction
+      //theTotalResult->Clear();
+      return theTotalResult;
+    }    
+  }
 
   // Next check for illegal track status
   //
diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc geant4.9.5.p01/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc
--- geant4.9.5.p01.orig/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc      2012-07-23 08:38:33.000000000 +0200
+++ geant4.9.5.p01/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc   2012-07-23 09:14:51.000000000 +0200
@@ -99,6 +99,9 @@ G4VParticleChange* G4WHadronElasticProce
                                  const G4Step& step)
 {
   aParticleChange.Initialize(track);
+  G4double weight = track.GetWeight();
+  aParticleChange.ProposeWeight(weight);
+
   G4double kineticEnergy = track.GetKineticEnergy();
   const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
   const G4ParticleDefinition* part = dynParticle->GetDefinition();
@@ -112,14 +115,14 @@ G4VParticleChange* G4WHadronElasticProce
       { return G4VDiscreteProcess::PostStepDoIt(track,step); }
 
   G4Material* material = track.GetMaterial();
-  G4Nucleus* targetNucleus = GetTargetNucleusPointer();
+  G4Nucleus* targNucleus = GetTargetNucleusPointer();
 
   // Select element
   G4Element* elm = 0;
   try
     {
       elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material, 
-                                                   *targetNucleus);
+                                                   *targNucleus);
     }
   catch(G4HadronicException & aR)
     {
@@ -138,8 +141,8 @@ G4VParticleChange* G4WHadronElasticProce
     {
       G4ExceptionDescription ed;
       ed << "Target element "<< elm->GetName()<<"  Z= " 
-        << targetNucleus->GetZ_asInt() << "  A= " 
-        << targetNucleus->GetA_asInt() << G4endl;
+        << targNucleus->GetZ_asInt() << "  A= " 
+        << targNucleus->GetA_asInt() << G4endl;
       DumpState(track,"ChooseHadronicInteraction",ed);
       ed << " No HadronicInteraction found out" << G4endl;
       G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had005", 
@@ -153,33 +156,38 @@ G4VParticleChange* G4WHadronElasticProce
 
   // Initialize the hadronic projectile from the track
   //  G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
-  G4HadProjectile thePro(track);
+  G4HadProjectile theProj(track);
   if(verboseLevel>1) {
     G4cout << "G4WHadronElasticProcess::PostStepDoIt for " 
           << part->GetParticleName()
           << " in " << material->GetName() 
-          << " Target Z= " << targetNucleus->GetZ_asInt() 
-          << " A= " << targetNucleus->GetA_asInt() << G4endl; 
+          << " Target Z= " << targNucleus->GetZ_asInt() 
+          << " A= " << targNucleus->GetA_asInt() << G4endl; 
   }
 
   G4HadFinalState* result = 0;
   try
     {
-      result = hadi->ApplyYourself( thePro, *targetNucleus);
+      result = hadi->ApplyYourself( theProj, *targNucleus);
     }
   catch(G4HadronicException aR)
     {
       G4ExceptionDescription ed;
       ed << "Call for " << hadi->GetModelName() << G4endl;
       ed << "Target element "<< elm->GetName()<<"  Z= " 
-        << targetNucleus->GetZ_asInt() 
-        << "  A= " << targetNucleus->GetA_asInt() << G4endl;
+        << targNucleus->GetZ_asInt() 
+        << "  A= " << targNucleus->GetA_asInt() << G4endl;
       DumpState(track,"ApplyYourself",ed);
       ed << " ApplyYourself failed" << G4endl;
       G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had006", 
                  FatalException, ed);
     }
 
+  // Check the result for catastrophic energy non-conservation
+  // cannot be applied because is not guranteed that recoil 
+  // nucleus is created
+  // result = CheckResult(theProj, targNucleus, result);
+
   // directions
   G4ThreeVector indir = track.GetMomentumDirection();
   G4double phi = CLHEP::twopi*G4UniformRand();
@@ -224,8 +232,9 @@ G4VParticleChange* G4WHadronElasticProce
 
   //G4cout << "Efinal= " << efinal << "  TrackStatus= " << status << G4endl;
 
-  // recoil
   aParticleChange.SetNumberOfSecondaries(0);
+
+  // recoil
   if(result->GetNumberOfSecondaries() > 0) {
     G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
 
@@ -238,7 +247,14 @@ G4VParticleChange* G4WHadronElasticProce
       pdir.rotateUz(indir);
       // G4cout << "recoil rotated " << pdir << G4endl;
       p->SetMomentumDirection(pdir);
-      aParticleChange.AddSecondary(p);
+
+      // in elastic scattering time and weight are not changed
+      G4Track* t = new G4Track(p, track.GetGlobalTime(), 
+                              track.GetPosition());
+      t->SetWeight(weight);
+      t->SetTouchableHandle(track.GetTouchableHandle());
+      aParticleChange.AddSecondary(t);
+
     } else {
       edep += p->GetKineticEnergy();
       delete p;
#	User	Rev	Content
1	davidlt	1.2	diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh geant4.9.5.p01/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh
2			--- geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh 2012-07-23 08:38:19.000000000 +0200
3			+++ geant4.9.5.p01/source/processes/hadronic/cross_sections/include/G4HadronCrossSections.hh 2012-07-23 09:14:23.000000000 +0200
4	davidlt	1.1	@@ -58,6 +58,8 @@
5
6			enum { TSIZE=41, NPARTS=35, NELAB=17, NCNLW=15, NFISS=21 };
7
8			+class G4Pow;
9			+
10			class G4HadronCrossSections
11			{
12			public:
13			@@ -100,6 +102,8 @@ class G4HadronCrossSections
14
15			static G4HadronCrossSections* theInstance;
16
17			+ G4Pow* g4pow;
18			+
19			G4double sigelastic;
20			G4double siginelastic;
21			G4ParticleDefinition* prevParticleDefinition;
22	davidlt	1.2	diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc geant4.9.5.p01/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc
23			--- geant4.9.5.p01.orig/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc 2012-07-23 08:38:19.000000000 +0200
24			+++ geant4.9.5.p01/source/processes/hadronic/cross_sections/src/G4HadronCrossSections.cc 2012-07-23 09:14:05.000000000 +0200
25	davidlt	1.1	@@ -55,7 +55,7 @@
26			#include "G4HadronCrossSections.hh"
27			#include "G4ios.hh"
28			#include "G4HadronicException.hh"
29			-
30			+#include "G4Pow.hh"
31
32			// Initialize static pointer for singleton instance
33			G4HadronCrossSections* G4HadronCrossSections::theInstance = 0;
34	davidlt	1.2	@@ -1217,7 +1217,10 @@ G4bool G4HadronCrossSections::correctIne
35	davidlt	1.1	G4HadronCrossSections::G4HadronCrossSections()
36			: prevParticleDefinition(0), prevElement(0), prevZZ(0), prevAA(0),
37			prevKineticEnergy(DBL_MAX), lastEkx(0.), lastEkxPower(0.), verboseLevel(0)
38			-{}
39			+{
40			+ siginelastic = sigelastic = 0.0;
41			+ g4pow = G4Pow::GetInstance();
42			+}
43
44			G4HadronCrossSections::~G4HadronCrossSections()
45			{}
46	davidlt	1.2	@@ -1243,7 +1246,12 @@ G4HadronCrossSections::GetInelasticCross
47	davidlt	1.1	if (particle->GetDefinition() != prevParticleDefinition \|\|
48			particle->GetKineticEnergy() != prevKineticEnergy \|\|
49			ZZ != prevZZ \|\| AA != prevAA)
50			- CalcScatteringCrossSections(particle, ZZ, AA);
51			+ { CalcScatteringCrossSections(particle, ZZ, AA); }
52			+ if (verboseLevel > 1) {
53			+ G4cout << " x-section inelastic(mb)= "
54			+ << siginelastic/millibarn << G4endl;
55			+ }
56			+
57			return siginelastic;
58			}
59
60	davidlt	1.2	@@ -1255,8 +1263,12 @@ G4HadronCrossSections::GetElasticCrossSe
61	davidlt	1.1	if (particle->GetDefinition() != prevParticleDefinition \|\|
62			particle->GetKineticEnergy() != prevKineticEnergy \|\|
63			ZZ != prevZZ \|\| AA != prevAA)
64			- CalcScatteringCrossSections(particle, ZZ, AA);
65			+ { CalcScatteringCrossSections(particle, ZZ, AA); }
66
67			+ if (verboseLevel > 1) {
68			+ G4cout << " x-section elastic(mb)= "
69			+ << sigelastic/millibarn << G4endl;
70			+ }
71			return sigelastic;
72			}
73
74	davidlt	1.2	@@ -1279,33 +1291,37 @@ G4HadronCrossSections::CalcScatteringCro
75	davidlt	1.1	G4double xspiel, xspiin;
76
77			G4int ipart = GetParticleCode(aParticle);
78			- G4double a = AA;
79			- G4double z = ZZ;
80			+ prevAA = AA;
81			+ prevZZ = ZZ;
82			+ prevParticleDefinition = aParticle->GetDefinition();
83			+ prevKineticEnergy = aParticle->GetKineticEnergy();
84
85			if (verboseLevel > 1) {
86			- G4cout << "G4HadronCrossSections: a=" << a << G4endl;
87			- G4cout << "G4HadronCrossSections: z=" << z << G4endl;
88			+ G4cout << "G4HadronCrossSections::ScatteringCrossSections: "
89			+ << aParticle->GetDefinition()->GetParticleName()
90			+ << " E(MeV)= " << aParticle->GetKineticEnergy()/MeV
91			+ << " Z= " << ZZ << " A= " << AA << G4endl;
92			}
93
94			-// Ions...
95			+ // Ions...
96
97			if (ipart >= 30 && ipart <= 32) {
98
99			G4double apart=0;
100			- if (ipart == 30) apart = std::pow(2., 1./3.);
101			- else if (ipart == 31) apart = std::pow(3., 1./3.);
102			- else if (ipart == 32) apart = std::pow(4., 1./3.);
103			+ if (ipart == 30) apart = g4pow->Z13(2);
104			+ else if (ipart == 31) apart = g4pow->Z13(3);
105			+ else if (ipart == 32) apart = g4pow->Z13(4);
106
107			- G4double term = apart + std::pow(a, 1./3.);
108			+ G4double term = apart + g4pow->Z13(AA);
109			sigin = 49.termterm;
110			- // Convert cross section from mb to default units
111			+ // Convert cross section from mb to default units
112			siginelastic = sigin*millibarn;
113			- if(aParticle->GetKineticEnergy() < 6*MeV) siginelastic = 0;
114			+ if(prevKineticEnergy < 6*MeV) siginelastic = 0;
115			sigelastic = 0.;
116			return;
117			}
118
119			- G4double ek = aParticle->GetKineticEnergy()/GeV;
120			+ G4double ek = prevKineticEnergy/GeV;
121
122			// Low energy neutrons...
123
124	davidlt	1.2	@@ -1349,7 +1365,7 @@ G4HadronCrossSections::CalcScatteringCro
125	davidlt	1.1	G4int ja2 = NCNLW - 1;
126			do {
127			G4int midBin = (ja1 + ja2)/2;
128			- if (a < cnlwat[midBin])
129			+ if (AA < cnlwat[midBin])
130			ja2 = midBin;
131			else
132			ja1 = midBin;
133	davidlt	1.2	@@ -1365,7 +1381,7 @@ G4HadronCrossSections::CalcScatteringCro
134	davidlt	1.1	dy = cnlwel[ja2][je1] - cnlwel[ja1][je1];
135			rca = dy/dnlwat;
136			b = cnlwel[ja1][je1] - rceelab[je1] - rcacnlwat[ja1];
137			- sigelastic = rceek + rcaa + b;
138			+ sigelastic = rceek + rcaAA + b;
139			// Inelastic cross section:
140			// E interpolation or extrapolation at JA1
141			dy = cnlwin[ja1][je2] - cnlwin[ja1][je1];
142	davidlt	1.2	@@ -1374,7 +1390,7 @@ G4HadronCrossSections::CalcScatteringCro
143	davidlt	1.1	dy = cnlwin[ja2][je1] - cnlwin[ja1][je1];
144			rca = dy/dnlwat;
145			b = cnlwin[ja1][je1] - rceelab[je1] - rcacnlwat[ja1];
146			- siginelastic = rceek + rcaa + b;
147			+ siginelastic = rceek + rcaAA + b;
148			// Convert cross sections from mb to default units
149			sigelastic = sigelastic*millibarn;
150			siginelastic = siginelastic*millibarn;
151	davidlt	1.2	@@ -1446,14 +1462,14 @@ G4HadronCrossSections::CalcScatteringCro
152	davidlt	1.1
153			// A-dependence from parameterization...
154
155			- if (a >= 1.5) {
156			+ if (AA >= 2) {
157
158			crel = 1.;
159			crin = 1.;
160
161	davidlt	1.2	G4int i = 3;
162	davidlt	1.1	- if (a < 50.) i = 2;
163			- if (a > 100.) i = 4;
164	davidlt	1.2	+ if (AA < 50) i = 2;
165			+ else if (AA > 100.) i = 4;
166			i = i - 1; // For array indexing
167	davidlt	1.1
168			// Protons and neutrons
169	davidlt	1.2	@@ -1474,7 +1490,7 @@ G4HadronCrossSections::CalcScatteringCro
170	davidlt	1.1	if (correctInelasticNearZero && je1 == 0 && i <= 3) {
171			G4double m0 = aParticle->GetMass()/GeV;
172			G4double T = std::sqrt(m0m0 + pp) - m0;
173			- G4double dx = std::sqrt(m0m0 + plab[1]plab[1]) - m0;
174			+ dx = std::sqrt(m0m0 + plab[1]plab[1]) - m0;
175			rc = dy/dx;
176			xsecin = rc*T + b;
177			}
178	davidlt	1.2	@@ -1490,7 +1506,7 @@ G4HadronCrossSections::CalcScatteringCro
179	davidlt	1.1	// Calculate correction factors (crel, crin) from values
180			// on Al, Cu, Pb. Note that data is only available for pions and protons.
181			G4double wgch = 0.5;
182			- if (a < 20.) wgch = 0.5 + 0.5*std::exp(-(a - 1.));
183			+ if (AA < 20) { wgch = 0.5 + 0.5*std::exp(-G4double(AA - 1)); }
184			sigel = wgchsigel + (1. - wgch)xsecel;
185			sigin = wgchsigin + (1. - wgch)xsecin;
186
187	davidlt	1.2	@@ -1517,8 +1533,8 @@ G4HadronCrossSections::CalcScatteringCro
188	davidlt	1.1	}
189
190			// Apply correction factors
191			- sigin = crin(sigin + sigel)std::pow(a, alph);
192			- sigel = crel0.36sigel*std::pow(a, 1.17);
193			+ sigin = crin(sigin + sigel)g4pow->powZ(AA, alph);
194			+ sigel = crel0.36sigel*g4pow->powZ(AA, 1.17);
195			sigel = sigel*partel[ipart1];
196			sigin = sigin*partin[ipart1];
197			}
198	davidlt	1.2	@@ -1591,7 +1607,7 @@ G4HadronCrossSections::GetFissionCrossSe
199	davidlt	1.1
200			G4double z43ba;
201			if (j == 4) {
202			- z43ba = std::pow(G4double(ZZ), 4./3.)/G4double(AA);
203			+ z43ba = ZZ*g4pow->Z13(ZZ)/G4double(AA);
204			z43ba = std::max(-67. + 38.7*z43ba, 0.);
205			} else {
206			z43ba = 1.;
207	davidlt	1.2	diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/management/src/G4HadronicProcess.cc geant4.9.5.p01/source/processes/hadronic/management/src/G4HadronicProcess.cc
208			--- geant4.9.5.p01.orig/source/processes/hadronic/management/src/G4HadronicProcess.cc 2012-07-23 08:38:33.000000000 +0200
209			+++ geant4.9.5.p01/source/processes/hadronic/management/src/G4HadronicProcess.cc 2012-07-23 09:15:11.000000000 +0200
210			@@ -213,11 +213,13 @@ G4HadronicProcess::PostStepDoIt(const G4
211			ed);
212			}
213
214			- if (GetElementCrossSection(aParticle, anElement, aMaterial) <= 0.0) {
215			- // No interaction
216			- //theTotalResult->Clear();
217			- return theTotalResult;
218			- }
219			+ if(aParticle->GetDefinition()->GetPDGCharge() != 0.0) {
220			+ if (GetElementCrossSection(aParticle, anElement, aMaterial) <= 0.0) {
221			+ // No interaction
222			+ //theTotalResult->Clear();
223			+ return theTotalResult;
224			+ }
225			+ }
226
227			// Next check for illegal track status
228			//
229			diff -arupN geant4.9.5.p01.orig/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc geant4.9.5.p01/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc
230			--- geant4.9.5.p01.orig/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc 2012-07-23 08:38:33.000000000 +0200
231			+++ geant4.9.5.p01/source/processes/hadronic/processes/src/G4WHadronElasticProcess.cc 2012-07-23 09:14:51.000000000 +0200
232			@@ -99,6 +99,9 @@ G4VParticleChange* G4WHadronElasticProce
233			const G4Step& step)
234			{
235			aParticleChange.Initialize(track);
236			+ G4double weight = track.GetWeight();
237			+ aParticleChange.ProposeWeight(weight);
238			+
239			G4double kineticEnergy = track.GetKineticEnergy();
240			const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
241			const G4ParticleDefinition* part = dynParticle->GetDefinition();
242			@@ -112,14 +115,14 @@ G4VParticleChange* G4WHadronElasticProce
243			{ return G4VDiscreteProcess::PostStepDoIt(track,step); }
244
245			G4Material* material = track.GetMaterial();
246			- G4Nucleus* targetNucleus = GetTargetNucleusPointer();
247			+ G4Nucleus* targNucleus = GetTargetNucleusPointer();
248
249			// Select element
250			G4Element* elm = 0;
251			try
252			{
253			elm = GetCrossSectionDataStore()->SampleZandA(dynParticle, material,
254			- *targetNucleus);
255			+ *targNucleus);
256			}
257			catch(G4HadronicException & aR)
258			{
259			@@ -138,8 +141,8 @@ G4VParticleChange* G4WHadronElasticProce
260			{
261			G4ExceptionDescription ed;
262			ed << "Target element "<< elm->GetName()<<" Z= "
263			- << targetNucleus->GetZ_asInt() << " A= "
264			- << targetNucleus->GetA_asInt() << G4endl;
265			+ << targNucleus->GetZ_asInt() << " A= "
266			+ << targNucleus->GetA_asInt() << G4endl;
267			DumpState(track,"ChooseHadronicInteraction",ed);
268			ed << " No HadronicInteraction found out" << G4endl;
269			G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had005",
270			@@ -153,33 +156,38 @@ G4VParticleChange* G4WHadronElasticProce
271
272			// Initialize the hadronic projectile from the track
273			// G4cout << "track " << track.GetDynamicParticle()->Get4Momentum()<<G4endl;
274			- G4HadProjectile thePro(track);
275			+ G4HadProjectile theProj(track);
276			if(verboseLevel>1) {
277			G4cout << "G4WHadronElasticProcess::PostStepDoIt for "
278			<< part->GetParticleName()
279			<< " in " << material->GetName()
280			- << " Target Z= " << targetNucleus->GetZ_asInt()
281			- << " A= " << targetNucleus->GetA_asInt() << G4endl;
282			+ << " Target Z= " << targNucleus->GetZ_asInt()
283			+ << " A= " << targNucleus->GetA_asInt() << G4endl;
284			}
285
286			G4HadFinalState* result = 0;
287			try
288			{
289			- result = hadi->ApplyYourself( thePro, *targetNucleus);
290			+ result = hadi->ApplyYourself( theProj, *targNucleus);
291			}
292			catch(G4HadronicException aR)
293			{
294			G4ExceptionDescription ed;
295			ed << "Call for " << hadi->GetModelName() << G4endl;
296			ed << "Target element "<< elm->GetName()<<" Z= "
297			- << targetNucleus->GetZ_asInt()
298			- << " A= " << targetNucleus->GetA_asInt() << G4endl;
299			+ << targNucleus->GetZ_asInt()
300			+ << " A= " << targNucleus->GetA_asInt() << G4endl;
301			DumpState(track,"ApplyYourself",ed);
302			ed << " ApplyYourself failed" << G4endl;
303			G4Exception("G4WHadronElasticProcess::PostStepDoIt", "had006",
304			FatalException, ed);
305			}
306
307			+ // Check the result for catastrophic energy non-conservation
308			+ // cannot be applied because is not guranteed that recoil
309			+ // nucleus is created
310			+ // result = CheckResult(theProj, targNucleus, result);
311			+
312			// directions
313			G4ThreeVector indir = track.GetMomentumDirection();
314			G4double phi = CLHEP::twopi*G4UniformRand();
315			@@ -224,8 +232,9 @@ G4VParticleChange* G4WHadronElasticProce
316
317			//G4cout << "Efinal= " << efinal << " TrackStatus= " << status << G4endl;
318
319			- // recoil
320			aParticleChange.SetNumberOfSecondaries(0);
321			+
322			+ // recoil
323			if(result->GetNumberOfSecondaries() > 0) {
324			G4DynamicParticle* p = result->GetSecondary(0)->GetParticle();
325
326			@@ -238,7 +247,14 @@ G4VParticleChange* G4WHadronElasticProce
327			pdir.rotateUz(indir);
328			// G4cout << "recoil rotated " << pdir << G4endl;
329			p->SetMomentumDirection(pdir);
330			- aParticleChange.AddSecondary(p);
331			+
332			+ // in elastic scattering time and weight are not changed
333			+ G4Track* t = new G4Track(p, track.GetGlobalTime(),
334			+ track.GetPosition());
335			+ t->SetWeight(weight);
336			+ t->SetTouchableHandle(track.GetTouchableHandle());
337			+ aParticleChange.AddSecondary(t);
338			+
339			} else {
340			edep += p->GetKineticEnergy();
341			delete p;