Index: trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiNeutronInelastic.cc
===================================================================
--- trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiNeutronInelastic.cc	(revision 1340)
+++ trunk/source/processes/hadronic/models/high_energy/src/G4HEAntiNeutronInelastic.cc	(revision 1347)
@@ -24,8 +24,6 @@
 // ********************************************************************
 //
-//
-// $Id: G4HEAntiNeutronInelastic.cc,v 1.15 2008/03/17 20:49:17 dennis Exp $
-// GEANT4 tag $Name: geant4-09-03-ref-09 $
-//
+// $Id: G4HEAntiNeutronInelastic.cc,v 1.17 2010/11/29 05:44:44 dennis Exp $
+// GEANT4 tag $Name: geant4-09-04-ref-00 $
 //
 
@@ -33,10 +31,9 @@
 #include "G4ios.hh"
 
-//
 // G4 Process: Gheisha High Energy Collision model.
 // This includes the high energy cascading model, the two-body-resonance model
-// and the low energy two-body model. Not included are the low energy stuff like
-// nuclear reactions, nuclear fission without any cascading and all processes for
-// particles at rest.  
+// and the low energy two-body model.  Not included is the low energy stuff
+// like nuclear reactions, nuclear fission without any cascading and all
+// processes for particles at rest.  
 // First work done by J.L.Chuma and F.W.Jones, TRIUMF, June 96.  
 // H. Fesefeldt, RWTH-Aachen, 23-October-1996
@@ -45,198 +42,177 @@
 #include "G4HEAntiNeutronInelastic.hh"
 
-G4HadFinalState *  G4HEAntiNeutronInelastic::
-ApplyYourself( const G4HadProjectile &aTrack, G4Nucleus &targetNucleus )
-  {
-    G4HEVector * pv = new G4HEVector[MAXPART];
-    const G4HadProjectile *aParticle = &aTrack;
-//    G4DynamicParticle *originalTarget = targetNucleus.ReturnTargetParticle();
-    const G4double atomicWeight = targetNucleus.GetN();
-    const G4double atomicNumber = targetNucleus.GetZ();
-    G4HEVector incidentParticle(aParticle);
-
-    G4int    incidentCode          = incidentParticle.getCode();
-    G4double incidentMass          = incidentParticle.getMass();
-    G4double incidentTotalEnergy   = incidentParticle.getEnergy();
-    G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
-    G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
-
-    if(incidentKineticEnergy < 1.)
-      { 
-        G4cout << "GHEAntiNeutronInelastic: incident energy < 1 GeV" << G4endl;;
-      }
-    if(verboseLevel > 1)
-      {
-        G4cout << "G4HEAntiNeutronInelastic::ApplyYourself" << G4endl;
-        G4cout << "incident particle " << incidentParticle.getName()
-             << "mass "              << incidentMass
-             << "kinetic energy "    << incidentKineticEnergy
-             << G4endl;
-        G4cout << "target material with (A,Z) = (" 
-             << atomicWeight << "," << atomicNumber << ")" << G4endl;
-      }
-
-    G4double inelasticity  = NuclearInelasticity(incidentKineticEnergy, 
-                                                 atomicWeight, atomicNumber);
-    if(verboseLevel > 1)
-        G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
+G4HadFinalState*
+G4HEAntiNeutronInelastic::ApplyYourself(const G4HadProjectile &aTrack,
+                                        G4Nucleus &targetNucleus)
+{
+  G4HEVector* pv = new G4HEVector[MAXPART];
+  const G4HadProjectile* aParticle = &aTrack;
+  const G4double atomicWeight = targetNucleus.GetN();
+  const G4double atomicNumber = targetNucleus.GetZ();
+  G4HEVector incidentParticle(aParticle);
+
+  G4int incidentCode = incidentParticle.getCode();
+  G4double incidentMass = incidentParticle.getMass();
+  G4double incidentTotalEnergy = incidentParticle.getEnergy();
+  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
+  G4double incidentKineticEnergy = incidentTotalEnergy - incidentMass;
+
+  if (incidentKineticEnergy < 1.)
+    G4cout << "GHEAntiNeutronInelastic: incident energy < 1 GeV" << G4endl;;
+
+  if (verboseLevel > 1) {
+    G4cout << "G4HEAntiNeutronInelastic::ApplyYourself" << G4endl;
+    G4cout << "incident particle " << incidentParticle.getName()
+           << "mass "              << incidentMass
+           << "kinetic energy "    << incidentKineticEnergy
+           << G4endl;
+    G4cout << "target material with (A,Z) = (" 
+           << atomicWeight << "," << atomicNumber << ")" << G4endl;
+  }
+
+  G4double inelasticity = NuclearInelasticity(incidentKineticEnergy, 
+                                              atomicWeight, atomicNumber);
+  if (verboseLevel > 1)
+    G4cout << "nuclear inelasticity = " << inelasticity << G4endl;
     
-    incidentKineticEnergy -= inelasticity;
+  incidentKineticEnergy -= inelasticity;
     
-    G4double excitationEnergyGNP = 0.;
-    G4double excitationEnergyDTA = 0.; 
-
-    G4double excitation    = NuclearExcitation(incidentKineticEnergy,
-                                               atomicWeight, atomicNumber,
-                                               excitationEnergyGNP,
-                                               excitationEnergyDTA);
-    if(verboseLevel > 1)
-      G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
+  G4double excitationEnergyGNP = 0.;
+  G4double excitationEnergyDTA = 0.; 
+
+  G4double excitation = NuclearExcitation(incidentKineticEnergy,
+                                          atomicWeight, atomicNumber,
+                                          excitationEnergyGNP,
+                                          excitationEnergyDTA);
+  if (verboseLevel > 1)
+    G4cout << "nuclear excitation = " << excitation << excitationEnergyGNP 
            << excitationEnergyDTA << G4endl;             
 
-
-    incidentKineticEnergy -= excitation;
-    incidentTotalEnergy    = incidentKineticEnergy + incidentMass;
-    incidentTotalMomentum  = std::sqrt( (incidentTotalEnergy-incidentMass)                    
-                                  *(incidentTotalEnergy+incidentMass));
-
-
-    G4HEVector targetParticle;
-    if(G4UniformRand() < atomicNumber/atomicWeight)
-      { 
-        targetParticle.setDefinition("Proton");
-      }
-    else
-      { 
-        targetParticle.setDefinition("Neutron");
-      }
-
-    G4double targetMass         = targetParticle.getMass();
-    G4double centerOfMassEnergy = std::sqrt( incidentMass*incidentMass + targetMass*targetMass
-                                       + 2.0*targetMass*incidentTotalEnergy);
-    G4double availableEnergy    = centerOfMassEnergy - targetMass - incidentMass;
-
-                                                                // this was the meaning of inElastic in the
-                                                                // original Gheisha stand-alone version. 
-//    G4bool   inElastic          = InElasticCrossSectionInFirstInt
-//                                    (availableEnergy, incidentCode, incidentTotalMomentum);  
-                                                                // by unknown reasons, it has been replaced
-                                                                // to the following code in Geant???
-    G4bool inElastic = true;
-//    if (G4UniformRand() < elasticCrossSection/totalCrossSection) inElastic = false;    
-
-    vecLength = 0;           
+  incidentKineticEnergy -= excitation;
+  incidentTotalEnergy = incidentKineticEnergy + incidentMass;
+  incidentTotalMomentum = std::sqrt( (incidentTotalEnergy-incidentMass)
+                                    *(incidentTotalEnergy+incidentMass));
+
+  G4HEVector targetParticle;
+  if (G4UniformRand() < atomicNumber/atomicWeight) { 
+    targetParticle.setDefinition("Proton");
+  } else { 
+    targetParticle.setDefinition("Neutron");
+  }
+
+  G4double targetMass = targetParticle.getMass();
+  G4double centerOfMassEnergy = std::sqrt(incidentMass*incidentMass
+                                        + targetMass*targetMass
+                                        + 2.0*targetMass*incidentTotalEnergy);
+  G4double availableEnergy = centerOfMassEnergy - targetMass - incidentMass;
+
+  G4bool inElastic = true;
+  vecLength = 0;           
         
-    if(verboseLevel > 1)
-      G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
+  if (verboseLevel > 1)
+    G4cout << "ApplyYourself: CallFirstIntInCascade for particle "
            << incidentCode << G4endl;
 
-    G4bool successful = false; 
+  G4bool successful = false; 
     
-    if(inElastic || (!inElastic && atomicWeight < 1.5))
-      { 
-        FirstIntInCasAntiNeutron(inElastic, availableEnergy, pv, vecLength,
-                                 incidentParticle, targetParticle, atomicWeight);
-
-        if(verboseLevel > 1)
-	   G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;  
-
-
-        if ((vecLength > 0) && (availableEnergy > 1.)) 
-                   StrangeParticlePairProduction( availableEnergy, centerOfMassEnergy,
-                                                  pv, vecLength,
-                                                  incidentParticle, targetParticle);
-            HighEnergyCascading( successful, pv, vecLength,
-                                 excitationEnergyGNP, excitationEnergyDTA,
-                                 incidentParticle, targetParticle,
-                                 atomicWeight, atomicNumber);
-        if (!successful)
-            HighEnergyClusterProduction( successful, pv, vecLength,
-                                         excitationEnergyGNP, excitationEnergyDTA,
-                                         incidentParticle, targetParticle,
-                                         atomicWeight, atomicNumber);
-        if (!successful) 
-            MediumEnergyCascading( successful, pv, vecLength, 
-                                   excitationEnergyGNP, excitationEnergyDTA, 
-                                   incidentParticle, targetParticle,
-                                   atomicWeight, atomicNumber);
-
-        if (!successful)
-            MediumEnergyClusterProduction( successful, pv, vecLength,
-                                           excitationEnergyGNP, excitationEnergyDTA,       
-                                           incidentParticle, targetParticle,
-                                           atomicWeight, atomicNumber);
-        if (!successful)
-            QuasiElasticScattering( successful, pv, vecLength,
-                                    excitationEnergyGNP, excitationEnergyDTA,
-                                    incidentParticle, targetParticle, 
-                                    atomicWeight, atomicNumber);
-      }
-    if (!successful)
-      { 
-            ElasticScattering( successful, pv, vecLength,
-                               incidentParticle,    
-                               atomicWeight, atomicNumber);
-      }
-
-    if (!successful)
-      { 
-        G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" << G4endl;
-      }
-      FillParticleChange(pv,  vecLength);
-      delete [] pv;
-      theParticleChange.SetStatusChange(stopAndKill);
-      return & theParticleChange;
-  }
+  FirstIntInCasAntiNeutron(inElastic, availableEnergy, pv, vecLength,
+                           incidentParticle, targetParticle, atomicWeight);
+
+  if (verboseLevel > 1)
+    G4cout << "ApplyYourself::StrangeParticlePairProduction" << G4endl;
+
+  if ((vecLength > 0) && (availableEnergy > 1.)) 
+    StrangeParticlePairProduction(availableEnergy, centerOfMassEnergy,
+                                  pv, vecLength,
+                                  incidentParticle, targetParticle);
+  HighEnergyCascading(successful, pv, vecLength,
+                      excitationEnergyGNP, excitationEnergyDTA,
+                      incidentParticle, targetParticle,
+                      atomicWeight, atomicNumber);
+  if (!successful)
+    HighEnergyClusterProduction(successful, pv, vecLength,
+                                excitationEnergyGNP, excitationEnergyDTA,
+                                incidentParticle, targetParticle,
+                                atomicWeight, atomicNumber);
+  if (!successful) 
+    MediumEnergyCascading(successful, pv, vecLength, 
+                          excitationEnergyGNP, excitationEnergyDTA, 
+                          incidentParticle, targetParticle,
+                          atomicWeight, atomicNumber);
+
+  if (!successful)
+    MediumEnergyClusterProduction(successful, pv, vecLength,
+                                  excitationEnergyGNP, excitationEnergyDTA,       
+                                  incidentParticle, targetParticle,
+                                  atomicWeight, atomicNumber);
+  if (!successful)
+    QuasiElasticScattering(successful, pv, vecLength,
+                           excitationEnergyGNP, excitationEnergyDTA,
+                           incidentParticle, targetParticle, 
+                           atomicWeight, atomicNumber);
+  if (!successful)
+    ElasticScattering(successful, pv, vecLength,
+                      incidentParticle,    
+                      atomicWeight, atomicNumber);
+
+  if (!successful)
+    G4cout << "GHEInelasticInteraction::ApplyYourself fails to produce final state particles" 
+           << G4endl;
+
+  FillParticleChange(pv,  vecLength);
+  delete [] pv;
+  theParticleChange.SetStatusChange(stopAndKill);
+  return &theParticleChange;
+}
+
 
 void
-  G4HEAntiNeutronInelastic::FirstIntInCasAntiNeutron( G4bool &inElastic,
-                                                      const G4double availableEnergy,
-                                                      G4HEVector pv[],
-                                                      G4int &vecLen,
-                                                      G4HEVector incidentParticle,
-                                                      G4HEVector targetParticle,
-                                                      const G4double atomicWeight)
-
-// AntiNeutron undergoes interaction with nucleon within a nucleus.  Check if it is
-// energetically possible to produce pions/kaons.  In not, assume nuclear excitation
-// occurs and input particle is degraded in energy. No other particles are produced.
+G4HEAntiNeutronInelastic::FirstIntInCasAntiNeutron(G4bool& inElastic,
+                                                   const G4double availableEnergy,
+                                                   G4HEVector pv[],
+                                                   G4int& vecLen,
+                                                   const G4HEVector& incidentParticle,
+                                                   const G4HEVector& targetParticle,
+                                                   const G4double atomicWeight)
+
+// AntiNeutron undergoes interaction with nucleon within a nucleus.  Check if
+// it is energetically possible to produce pions/kaons.  If not, assume
+// nuclear excitation occurs and input particle is degraded in energy.  No
+// other particles are produced.
 // If reaction is possible, find the correct number of pions/protons/neutrons
 // produced using an interpolation to multiplicity data.  Replace some pions or
 // protons/neutrons by kaons or strange baryons according to the average
 // multiplicity per inelastic reaction.
-
- {
-   static const G4double expxu =  std::log(MAXFLOAT); // upper bound for arg. of exp
-   static const G4double expxl = -expxu;         // lower bound for arg. of exp
-
-   static const G4double protb = 0.7;
-   static const G4double neutb = 0.7;
-   static const G4double     c = 1.25;
-
-   static const G4int   numMul   = 1200;
-   static const G4int   numMulAn = 400;
-   static const G4int   numSec   = 60;
-
-   G4int              neutronCode = Neutron.getCode();
-   G4int              protonCode  = Proton.getCode();
-
-   G4int               targetCode = targetParticle.getCode();
-//   G4double          incidentMass = incidentParticle.getMass();
-//   G4double        incidentEnergy = incidentParticle.getEnergy();
-   G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
-
-   static G4bool first = true;
-   static G4double protmul[numMul],  protnorm[numSec];   // proton constants
-   static G4double protmulAn[numMulAn],protnormAn[numSec]; 
-   static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
-   static G4double neutmulAn[numMulAn],neutnormAn[numSec];
-
-                              //  misc. local variables
-                              //  np = number of pi+,  nm = number of pi-,  nz = number of pi0
-
-   G4int i, counter, nt, np, nm, nz;
+{
+  static const G4double expxu = std::log(MAXFLOAT); // upper bound for arg. of exp
+  static const G4double expxl = -expxu;             // lower bound for arg. of exp
+
+  static const G4double protb = 0.7;
+  static const G4double neutb = 0.7;
+  static const G4double     c = 1.25;
+
+  static const G4int numMul   = 1200;
+  static const G4int numMulAn = 400;
+  static const G4int numSec   = 60;
+
+  G4int neutronCode = Neutron.getCode();
+  G4int protonCode = Proton.getCode();
+
+  G4int targetCode = targetParticle.getCode();
+  G4double incidentTotalMomentum = incidentParticle.getTotalMomentum();
+
+  static G4bool first = true;
+  static G4double protmul[numMul],  protnorm[numSec];   // proton constants
+  static G4double protmulAn[numMulAn],protnormAn[numSec]; 
+  static G4double neutmul[numMul],  neutnorm[numSec];   // neutron constants
+  static G4double neutmulAn[numMulAn],neutnormAn[numSec];
+
+  //  misc. local variables
+  //  np = number of pi+,  nm = number of pi-,  nz = number of pi0
+
+  G4int i, counter, nt, np, nm, nz;
 
    if( first ) 
-     {                         // compute normalization constants, this will only be done once
+     {                 // compute normalization constants, this will only be done once
        first = false;
        for( i=0; i<numMul  ; i++ ) protmul[i]  = 0.0;
