Index: trunk/source/processes/hadronic/models/cascade/cascade/src/G4NucleiModel.cc
===================================================================
--- trunk/source/processes/hadronic/models/cascade/cascade/src/G4NucleiModel.cc	(revision 1337)
+++ trunk/source/processes/hadronic/models/cascade/cascade/src/G4NucleiModel.cc	(revision 1340)
@@ -23,7 +23,6 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
-//
-// $Id: G4NucleiModel.cc,v 1.48.2.1 2010/06/25 09:44:56 gunter Exp $
-// Geant4 tag: $Name: geant4-09-04-beta-01 $
+// $Id: G4NucleiModel.cc,v 1.87 2010/09/24 06:26:06 mkelsey Exp $
+// Geant4 tag: $Name: geant4-09-03-ref-09 $
 //
 // 20100112  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
@@ -41,8 +40,34 @@
 //		calculations.  use G4Cascade*Channel for total xsec in pi-N
 //		N-N channels.  Move absorptionCrossSection() from SpecialFunc.
-
-//#define CHC_CHECK
+// 20100610  M. Kelsey -- Replace another random-angle code block; add some
+//		diagnostic output for partner-list production.
+// 20100617  M. Kelsey -- Replace preprocessor flag CHC_CHECK with
+//		G4CASCADE_DEBUG_CHARGE
+// 20100620  M. Kelsey -- Improve error message on empty partners list, add
+//		four-momentum checking after EPCollider
+// 20100621  M. Kelsey -- In boundaryTransition() account for momentum transfer
+//		to secondary by boosting into recoil nucleus "rest" frame.
+//		Don't need temporaries to create dummy "partners" for list.
+// 20100622  M. Kelsey -- Restore use of "bindingEnergy()" function name, which
+//		is now a wrapper for G4NucleiProperties::GetBindingEnergy().
+// 20100623  M. Kelsey -- Eliminate some temporaries terminating partner-list,
+//		discard recoil boost for now. Initialize all data
+//		members in ctors.  Allow generateModel() to be called
+//		mutliple times, by clearing vectors each time through;
+//		avoid extra work by returning if A and Z are same as
+//		before.
+// 20100628  M. Kelsey -- Two momentum-recoil bugs; don't subtract energies!
+// 20100715  M. Kelsey -- Make G4InuclNuclei in generateModel(), use for
+//		balance checking (only verbose>2) in generateParticleFate().
+// 20100721  M. Kelsey -- Use new G4CASCADE_CHECK_ECONS for conservation checks
+// 20100723  M. Kelsey -- Move G4CollisionOutput buffer to .hh for reuse
+// 20100726  M. Kelsey -- Preallocate arrays with number_of_zones dimension.
+// 20100902  M. Kelsey -- Remove resize(3) directives from qdeutron/acsecs
+// 20100907  M. Kelsey -- Limit interaction targets based on current nucleon
+//		counts (e.g., no pp if protonNumberCurrent < 2!)
+// 20100914  M. Kelsey -- Migrate to integer A and Z
 
 #include "G4NucleiModel.hh"
+#include "G4CascadeCheckBalance.hh"
 #include "G4CascadeInterpolator.hh"
 #include "G4CascadeNNChannel.hh"
@@ -63,5 +88,4 @@
 #include "G4LorentzConvertor.hh"
 #include "G4Neutron.hh"
-#include "G4NucleiProperties.hh"
 #include "G4Proton.hh"
 
@@ -69,14 +93,31 @@
 using namespace G4InuclSpecialFunctions;
 
-
 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
 
-G4NucleiModel::G4NucleiModel() : verboseLevel(0) {
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4NucleiModel::G4NucleiModel" << G4endl;
-  }
-}
-
-G4NucleiModel::G4NucleiModel(G4InuclNuclei* nuclei) : verboseLevel(0) {
+G4NucleiModel::G4NucleiModel()
+  : verboseLevel(0), A(0), Z(0), theNucleus(0),
+    neutronNumber(0), protonNumber(0),
+    neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
+    current_nucl2(0) {}
+
+G4NucleiModel::G4NucleiModel(G4int a, G4int z)
+  : verboseLevel(0), A(0), Z(0), theNucleus(0),
+    neutronNumber(0), protonNumber(0),
+    neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
+    current_nucl2(0) {
+  generateModel(a,z);
+}
+
+G4NucleiModel::G4NucleiModel(G4InuclNuclei* nuclei)
+  : verboseLevel(0), A(0), Z(0), theNucleus(0),
+    neutronNumber(0), protonNumber(0),
+    neutronNumberCurrent(0), protonNumberCurrent(0), current_nucl1(0),
+    current_nucl2(0) {
+  generateModel(nuclei);
+}
+
+
+void 
+G4NucleiModel::generateModel(G4InuclNuclei* nuclei) {
   generateModel(nuclei->getA(), nuclei->getZ());
 }
@@ -84,7 +125,8 @@
 
 void 
-G4NucleiModel::generateModel(G4double a, G4double z) {
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4NucleiModel::generateModel" << G4endl;
+G4NucleiModel::generateModel(G4int a, G4int z) {
+  if (verboseLevel) {
+    G4cout << " >>> G4NucleiModel::generateModel A " << a << " Z " << z
+	   << G4endl;
   }
 
@@ -95,6 +137,6 @@
   const G4double pion_vp = 0.007; // in GeV
   const G4double pion_vp_small = 0.007; 
-  const G4double radForSmall = 8.0; // fermi
-  const G4double piTimes4thirds = 4.189; // 4 Pi/3
+  const G4double radForSmall = 8.0;	// Units 0.f fm?  R_p = 0.8768 fm
+  const G4double piTimes4thirds = 4.189; // 4 Pi/3	FIXME!
   const G4double mproton = G4Proton::Definition()->GetPDGMass() / GeV;
   const G4double mneutron = G4Neutron::Definition()->GetPDGMass() / GeV;
@@ -102,17 +144,34 @@
   const G4double alfa6[6] = { 0.9, 0.6, 0.4, 0.2, 0.1, 0.05 };
 
+  // If model already built, just return; otherwise intialize everything
+  if (a == A && z == Z) {
+    if (verboseLevel > 1)
+      G4cout << " model already generated for A=" << a << ", Z=" << z << G4endl;
+
+    reset();		// Zeros out neutron/proton evaporates
+    return;
+  }
+
   A = a;
   Z = z;
+  delete theNucleus;
+  theNucleus = new G4InuclNuclei(A,Z);		// For conservation checking
+
   neutronNumber = a - z;
   protonNumber = z;
-  neutronNumberCurrent = neutronNumber;
-  protonNumberCurrent = protonNumber;
-
-// Set binding energies
-//  G4double dm = bindingEnergy(a, z);
-  G4double dm = G4NucleiProperties::GetBindingEnergy(G4lrint(a), G4lrint(z));
-
-  binding_energies.push_back(0.001 * std::fabs(G4NucleiProperties::GetBindingEnergy(G4lrint(a-1), G4lrint(z-1)) - dm)); // for P
-  binding_energies.push_back(0.001 * std::fabs(G4NucleiProperties::GetBindingEnergy(G4lrint(a-1), G4lrint(z)) - dm)); // for N
+  reset();
+
+  // Clear all parameters arrays for reloading
+  binding_energies.clear();
+  nucleon_densities.clear();
+  zone_potentials.clear();
+  fermi_momenta.clear();
+  zone_radii.clear();
+
+  // Set binding energies
+  G4double dm = bindingEnergy(a,z);
+
+  binding_energies.push_back(0.001 * std::fabs(bindingEnergy(a-1,z-1)-dm)); // for P
+  binding_energies.push_back(0.001 * std::fabs(bindingEnergy(a-1,z)-dm));   // for N
 
   G4double CU = cuu*G4cbrt(a); // half-density radius * 2.8197
@@ -121,47 +180,49 @@
   G4double CU2 = 0.0; 
 
-  if (a > 4.5) {
-    std::vector<G4double> ur;
+  // This will be used to pre-allocate lots of arrays below
+  number_of_zones = (a < 5) ? 1 : (a < 100) ? 3 : 6;
+
+  // Buffers used for zone-by-zone calculations
+  G4double ur[7];
+  G4double v[6];
+  G4double v1[6];
+
+  // These are passed into nested vectors, so can't be made C arrays
+  std::vector<G4double> rod(number_of_zones);
+  std::vector<G4double> pf(number_of_zones);
+  std::vector<G4double> vz(number_of_zones);
+
+  if (a > 4) {
     G4int icase = 0;
 
-    if (a > 99.5) {
-      number_of_zones = 6;
-      ur.push_back(-D1);
-
+    if (a > 99) {
+      ur[0] = -D1;
       for (G4int i = 0; i < number_of_zones; i++) {
         G4double y = std::log((1.0 + D) / alfa6[i] - 1.0);
         zone_radii.push_back(CU + AU * y);
-        ur.push_back(y);
+        ur[i+1] = y;
       }
-
-    } else if (a > 11.5) {
-      number_of_zones = 3;
-      ur.push_back(-D1);
-
+    } else if (a > 11) {
+      ur[0] = -D1;
       for (G4int i = 0; i < number_of_zones; i++) {
 	G4double y = std::log((1.0 + D)/alfa3[i] - 1.0);
 	zone_radii.push_back(CU + AU * y);
-	ur.push_back(y);
+	ur[i+1] = y;
       }
-
     } else {
-      number_of_zones = 3;
       icase = 1;
-      ur.push_back(0.0);
  
       G4double CU1 = CU * CU;
       CU2 = std::sqrt(CU1 * (1.0 - 1.0 / a) + 6.4);
 
+      ur[0] = 0.0;
       for (G4int i = 0; i < number_of_zones; i++) {
 	G4double y = std::sqrt(-std::log(alfa3[i]));
 	zone_radii.push_back(CU2 * y);
-	ur.push_back(y);
+	ur[i+1] = y;
       }
     }
 
     G4double tot_vol = 0.0;
-    std::vector<G4double> v;
-    std::vector<G4double> v1;
-
     G4int i(0);
     for (i = 0; i < number_of_zones; i++) {
@@ -169,11 +230,10 @@
 
       if (icase == 0) {
-	v0 = volNumInt(ur[i], ur[i + 1], CU, D1);
-
+	v0 = volNumInt(ur[i], ur[i+1], CU, D1);
       } else {
-	v0 = volNumInt1(ur[i], ur[i + 1], CU2);
-      };
+	v0 = volNumInt1(ur[i], ur[i+1], CU2);
+      }
  
-      v.push_back(v0);
+      v[i] = v0;
       tot_vol += v0;
 
@@ -181,12 +241,12 @@
       if (i > 0) v0 -= zone_radii[i-1]*zone_radii[i-1]*zone_radii[i-1];
 
-      v1.push_back(v0);
+      v1[i] = v0;
     }
 
     // Protons
-    G4double dd0 = z/tot_vol/piTimes4thirds;
-    std::vector<G4double> rod;
-    std::vector<G4double> pf;
-    std::vector<G4double> vz;
+    G4double dd0 = G4double(z)/tot_vol/piTimes4thirds;
+    rod.clear();
+    pf.clear();
+    vz.clear();
 
     for (i = 0; i < number_of_zones; i++) {
@@ -203,5 +263,5 @@
 
     // Neutrons
-    dd0 = (a - z)/tot_vol/piTimes4thirds;
+    dd0 = G4double(a-z)/tot_vol/piTimes4thirds;
     rod.clear();
     pf.clear();
@@ -214,5 +274,5 @@
       pf.push_back(pff);
       vz.push_back(0.5 * pff * pff / mneutron + binding_energies[1]);
-    };
+    }
 
     nucleon_densities.push_back(rod);
@@ -221,17 +281,16 @@
 
     // pion stuff (primitive)
-    std::vector<G4double> vp(number_of_zones, pion_vp);
+    const std::vector<G4double> vp(number_of_zones, pion_vp);
     zone_potentials.push_back(vp);
 
     // kaon potential (primitive)
-    std::vector<G4double> kp(number_of_zones, -0.015);
+    const std::vector<G4double> kp(number_of_zones, -0.015);
     zone_potentials.push_back(kp);
 
     // hyperon potential (primitive)
-    std::vector<G4double> hp(number_of_zones, 0.03);
+    const std::vector<G4double> hp(number_of_zones, 0.03);
     zone_potentials.push_back(hp);
 
   } else { // a < 5
-    number_of_zones = 1;
     G4double smallRad = radForSmall;
     if (a == 4) smallRad *= 0.7;
@@ -240,7 +299,4 @@
 
     // proton
-    std::vector<G4double> rod;
-    std::vector<G4double> pf;
-    std::vector<G4double> vz;
     for (G4int i = 0; i < number_of_zones; i++) {
       G4double rd = vol*z;
@@ -272,58 +328,17 @@
 
     // pion (primitive)
-    std::vector<G4double> vp(number_of_zones, pion_vp_small);
+    const std::vector<G4double> vp(number_of_zones, pion_vp_small);
     zone_potentials.push_back(vp);
   
     // kaon potential (primitive)
-    std::vector<G4double> kp(number_of_zones, -0.015);
+    const std::vector<G4double> kp(number_of_zones, -0.015);
     zone_potentials.push_back(kp);
 
     // hyperon potential (primitive)
-    std::vector<G4double> hp(number_of_zones, 0.03);
+    const std::vector<G4double> hp(number_of_zones, 0.03);
     zone_potentials.push_back(hp);
   }
 
   nuclei_radius = zone_radii[zone_radii.size() - 1];
-
-  /*
-  // Print nuclear radii and densities
-  G4cout << " For A = " << a << " zone radii = ";
-  for (G4int i = 0; i < number_of_zones; i++) G4cout << zone_radii[i] << "  ";
-  G4cout << "  " << G4endl;
-
-  G4cout << " proton densities: ";
-  for (G4int i = 0; i < number_of_zones; i++)
-     G4cout << nucleon_densities[0][i] << "  ";
-  G4cout << G4endl;
-
-  G4cout << " neutron densities: ";
-  for (G4int i = 0; i < number_of_zones; i++)
-     G4cout << nucleon_densities[1][i] << "  ";
-  G4cout << G4endl;
-
-  G4cout << " protons per shell " ;
-  G4double rinner = 0.0;
-  G4double router = 0.0;
-  G4double shellVolume = 0.0;
-  for (G4int i = 0; i < number_of_zones; i++) {  // loop over zones
-    router = zone_radii[i];
-    shellVolume = piTimes4thirds*(router*router*router - rinner*rinner*rinner);
-    G4cout << G4lrint(shellVolume*nucleon_densities[0][i]) << "  ";
-    rinner = router;
-  }
-  G4cout << G4endl;
-
-  G4cout << " neutrons per shell " ;
-  rinner = 0.0;
-  router = 0.0;
-  shellVolume = 0.0;
-  for (G4int i = 0; i < number_of_zones; i++) {  // loop over zones
-    router = zone_radii[i];
-    shellVolume = piTimes4thirds*(router*router*router - rinner*rinner*rinner);
-    G4cout << G4lrint(shellVolume*nucleon_densities[1][i]) << "  ";
-    rinner = router;
-  }
-  G4cout << G4endl;
-  */
 }
 
@@ -345,6 +360,5 @@
 G4NucleiModel::volNumInt(G4double r1, G4double r2, 
 			 G4double, G4double d1) const {
-
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::volNumInt" << G4endl;
   }
@@ -373,5 +387,5 @@
       r += dr1; 
       fi += r * (r + d2) / (1.0 + std::exp(r));
-    };
+    }
 
     fun = 0.5 * fun1 + fi * dr;
@@ -398,5 +412,5 @@
 G4NucleiModel::volNumInt1(G4double r1, G4double r2, 
 			  G4double cu2) const {
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::volNumInt1" << G4endl;
   }
@@ -447,17 +461,15 @@
 
 void G4NucleiModel::printModel() const {
-
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::printModel" << G4endl;
   }
 
   G4cout << " nuclei model for A " << A << " Z " << Z << G4endl
-	 << " proton binding energy " << binding_energies[0] << 
-    " neutron binding energy " << binding_energies[1] << G4endl
-	 << " Nculei radius " << nuclei_radius << " number of zones " <<
-    number_of_zones << G4endl;
+	 << " proton binding energy " << binding_energies[0]
+	 << " neutron binding energy " << binding_energies[1] << G4endl
+	 << " Nuclei radius " << nuclei_radius << " number of zones "
+	 << number_of_zones << G4endl;
 
   for (G4int i = 0; i < number_of_zones; i++)
-
     G4cout << " zone " << i+1 << " radius " << zone_radii[i] << G4endl
 	   << " protons: density " << getDensity(1,i) << " PF " << 
@@ -480,5 +492,5 @@
 G4InuclElementaryParticle 
 G4NucleiModel::generateNucleon(G4int type, G4int zone) const {
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::generateNucleon" << G4endl;
   }
@@ -492,6 +504,5 @@
 G4NucleiModel::generateQuasiDeutron(G4int type1, G4int type2,
 				    G4int zone) const {
-
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::generateQuasiDeutron" << G4endl;
   }
@@ -518,9 +529,9 @@
 void
 G4NucleiModel::generateInteractionPartners(G4CascadParticle& cparticle) {
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::generateInteractionPartners" << G4endl;
   }
 
-  const G4double pi4by3 = 4.1887903; // 4 Pi / 3
+  const G4double pi4by3 = 4.1887903; // 4 Pi / 3	FIXME!
   const G4double small = 1.0e-10;
   const G4double huge_num = 50.0;
@@ -556,216 +567,227 @@
     r_in = zone_radii[zone - 1];
     r_out = zone_radii[zone];
-  };  
+  }  
 
   G4double path = cparticle.getPathToTheNextZone(r_in, r_out);
 
-  if (verboseLevel > 2){
+  if (verboseLevel > 2) {
     G4cout << " r_in " << r_in << " r_out " << r_out << " path " << path << G4endl;
   }
 
-  if (path < -small) { // something wrong
+  if (path < -small) { 			// something wrong
+    G4cerr << " generateInteractionPartners-> negative path length" << G4endl;
     return;
-
-  } else if (std::fabs(path) < small) { // just on the boundary
-    path = 0.0; 
-
-    G4InuclElementaryParticle particle;	// Dummy -- no type or momentum
-    thePartners.push_back(partner(particle, path));
-
-  } else { // normal case  
-    G4LorentzConvertor dummy_convertor;
-    dummy_convertor.setBullet(pmom, pmass);
-  
-    for (G4int ip = 1; ip < 3; ip++) { 
-      G4InuclElementaryParticle particle = generateNucleon(ip, zone);
-      dummy_convertor.setTarget(particle.getMomentum(), particle.getMass());
+  }
+
+  if (std::fabs(path) < small) { 	// just on the boundary
+    if (verboseLevel > 3)
+      G4cout << " generateInteractionPartners-> zero path" << G4endl;
+
+    thePartners.push_back(partner());	// Dummy list terminator with zero path
+    return;
+  }
+
+  G4LorentzConvertor dummy_convertor;
+  dummy_convertor.setBullet(pmom, pmass);
+
+  for (G4int ip = 1; ip < 3; ip++) {
+    // Only process nucleons which remain active in target
+    if (ip==1 && protonNumberCurrent < 1) continue;
+    if (ip==2 && neutronNumberCurrent < 1) continue;
+
+    // All nucleons are assumed to be at rest when colliding
+    G4InuclElementaryParticle particle = generateNucleon(ip, zone);
+    dummy_convertor.setTarget(particle.getMomentum(), particle.getMass());
+    G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
+    
+    G4double csec = totalCrossSection(ekin, ptype * ip);
+    
+    if(verboseLevel > 2) {
+      G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl;
+    }
+    
+    G4double dens = nucleon_densities[ip - 1][zone];
+    G4double rat = getRatio(ip);
+    G4double pw = -path * dens * csec * rat;
+    
+    if (pw < -huge_num) pw = -huge_num;
+    pw = 1.0 - std::exp(pw);
+    
+    if (verboseLevel > 2) {
+      G4cout << " pw " << pw << " rat " << rat << G4endl;
+    }
+    
+    G4double spath = path;
+    
+    if (inuclRndm() < pw) {
+      spath = -1.0 / dens / csec / rat * std::log(1.0 - pw * inuclRndm());
+      if (cparticle.young(young_cut, spath)) spath = path;
+      
+      if (verboseLevel > 2) {
+	G4cout << " ip " << ip << " spath " << spath << G4endl;
+      }
+    }
+
+    if (spath < path) {
+      if (verboseLevel > 3) {
+	G4cout << " adding partner[" << thePartners.size() << "]: ";
+	particle.printParticle();
+      }
+      thePartners.push_back(partner(particle, spath));
+    }
+  }  
+  
+  if (verboseLevel > 2) {
+    G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
+  }
+  
+  if (cparticle.getParticle().pion()) { // absorption possible
+    if (verboseLevel > 2) {
+      G4cout << " trying quasi-deuterons with bullet: ";
+      cparticle.getParticle().printParticle();
+    }
+
+    // Initialize buffers for quasi-deuteron results
+    qdeutrons.clear();
+    acsecs.clear();
+
+    G4double tot_abs_csec = 0.0;
+    G4double abs_sec;
+    G4double vol = zone_radii[zone]*zone_radii[zone]*zone_radii[zone];
+    
+    if (zone > 0) vol -= zone_radii[zone-1]*zone_radii[zone-1]*zone_radii[zone-1];
+    vol *= pi4by3; 
+    
+    G4double rat  = getRatio(1); 
+    G4double rat1 = getRatio(2); 
+
+    // FIXME:  Shouldn't be creating zero-cross-section states!
+
+    // Proton-proton state interacts with pi-, pi0 only
+    G4InuclElementaryParticle ppd = generateQuasiDeutron(1, 1, zone);
+    
+    if (protonNumberCurrent < 2 || !(ptype == pi0 || ptype == pim)) {
+      abs_sec = 0.0;
+    } else {
+      dummy_convertor.setTarget(ppd.getMomentum(), ppd.getMass());
+      
       G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
-
-      // Total cross section converted from mb to fm**2
-      G4double csec = totalCrossSection(ekin, ptype * ip);
-
-      if(verboseLevel > 2){
-	G4cout << " ip " << ip << " ekin " << ekin << " csec " << csec << G4endl;
+      
+      if (verboseLevel > 2) {
+	G4cout << " ptype=" << ptype << " using pp target" << G4endl;
+	ppd.printParticle();
       }
-
-      G4double dens = nucleon_densities[ip - 1][zone];
-      G4double rat = getRatio(ip);
-      G4double pw = -path * dens * csec * rat;
-
-      if (pw < -huge_num) pw = -huge_num;
-      pw = 1.0 - std::exp(pw);
-
-      if (verboseLevel > 2){
-	G4cout << " pw " << pw << " rat " << rat << G4endl;
-      }
-
-      G4double spath = path;
-
-      if (inuclRndm() < pw) {
-	spath = -1.0 / dens / csec / rat * std::log(1.0 - pw * inuclRndm());
-	if (cparticle.young(young_cut, spath)) spath = path;
-
-	if (verboseLevel > 2){
-	  G4cout << " ip " << ip << " spath " << spath << G4endl;
-	}
-
-      };
-      if (spath < path) thePartners.push_back(partner(particle, spath));
-    };  
-
-    if (verboseLevel > 2){
-      G4cout << " after nucleons " << thePartners.size() << " path " << path << G4endl;
-    }
-
-    if (cparticle.getParticle().pion()) { // absorption possible
-      if (verboseLevel > 2) {
-	G4cout << " trying quasi-deuterons with bullet: ";
-	cparticle.getParticle().printParticle();
-      }
-
-      std::vector<G4InuclElementaryParticle> qdeutrons(3);
-      std::vector<G4double> acsecs(3);
-
-      G4double tot_abs_csec = 0.0;
-      G4double abs_sec;
-      G4double vol = zone_radii[zone]*zone_radii[zone]*zone_radii[zone];
-
-      if (zone > 0) vol -= zone_radii[zone-1]*zone_radii[zone-1]*zone_radii[zone-1];
-      vol *= pi4by3; 
-
-      G4double rat  = getRatio(1); 
-      G4double rat1 = getRatio(2); 
-
-      G4InuclElementaryParticle ppd = generateQuasiDeutron(1, 1, zone);
-
-      if (ptype == 7 || ptype == 5) {
-	dummy_convertor.setTarget(ppd.getMomentum(), ppd.getMass());
-
-	G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
-
-	if (verboseLevel > 2) {
-	  G4cout << " ptype=" << ptype << " using pp target" << G4endl;
-	  ppd.printParticle();
-	}
-
-	abs_sec = absorptionCrossSection(ekin, ptype);
-	abs_sec *= nucleon_densities[0][zone] * nucleon_densities[0][zone]*
-	  rat * rat * vol; 
-
-      } else {
-	abs_sec = 0.0;
-      }; 
-
-      // abs_sec = 0.0;
-      tot_abs_csec += abs_sec;
-      acsecs.push_back(abs_sec);
-      qdeutrons.push_back(ppd);
-
-      G4InuclElementaryParticle npd = generateQuasiDeutron(1, 2, zone);
-
+      
+      abs_sec = absorptionCrossSection(ekin, ptype);
+      abs_sec *= nucleon_densities[0][zone] * nucleon_densities[0][zone]*
+	rat * rat * vol; 
+    }
+    
+    tot_abs_csec += abs_sec;
+    acsecs.push_back(abs_sec);
+    qdeutrons.push_back(ppd);
+    
+    // Proton-neutron state interacts with any pion type
+    G4InuclElementaryParticle npd = generateQuasiDeutron(1, 2, zone);
+
+    if (protonNumberCurrent < 1 || neutronNumberCurrent < 1) {
+      abs_sec = 0.0;
+    } else {
       dummy_convertor.setTarget(npd.getMomentum(), npd.getMass());
-
+      
       G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
-
+      
       if (verboseLevel > 2) {
 	G4cout << " using np target" << G4endl;
 	npd.printParticle();
       }
-
+      
       abs_sec = absorptionCrossSection(ekin, ptype); 
       abs_sec *= pn_spec * nucleon_densities[0][zone] * nucleon_densities[1][zone] *
-	rat * rat1 * vol; 
-      tot_abs_csec += abs_sec;
-      acsecs.push_back(abs_sec);
-      qdeutrons.push_back(npd);
-
-      G4InuclElementaryParticle nnd = generateQuasiDeutron(2, 2, zone);
-
-      if (ptype == 7 || ptype == 3) {
-	dummy_convertor.setTarget(nnd.getMomentum(), nnd.getMass());
-
-	G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
-
-	if (verboseLevel > 2) {
-	  G4cout << " ptype=" << ptype << " using nn target" << G4endl;
-	  nnd.printParticle();
+	rat * rat1 * vol;
+    }
+
+    tot_abs_csec += abs_sec;
+    acsecs.push_back(abs_sec);
+    qdeutrons.push_back(npd);
+
+    // Neutron-neutron state interacts with pi+, pi0 only
+    G4InuclElementaryParticle nnd = generateQuasiDeutron(2, 2, zone);
+    
+    if (neutronNumberCurrent < 2 || !(ptype == pi0 || ptype == pip)) {
+      abs_sec = 0.0;
+    } else {
+      dummy_convertor.setTarget(nnd.getMomentum(), nnd.getMass());
+      
+      G4double ekin = dummy_convertor.getKinEnergyInTheTRS();
+      
+      if (verboseLevel > 2) {
+	G4cout << " ptype=" << ptype << " using nn target" << G4endl;
+	nnd.printParticle();
+      }
+      
+      abs_sec = absorptionCrossSection(ekin, ptype); 
+      abs_sec *= nucleon_densities[1][zone] * nucleon_densities[1][zone] *
+	rat1 * rat1 * vol; 
+    } 
+
+    tot_abs_csec += abs_sec;
+    acsecs.push_back(abs_sec);
+    qdeutrons.push_back(nnd);
+
+    // Select quasideuteron interaction from non-zero cross-section choices
+    if (verboseLevel > 2){
+      G4cout << " rod1 " << acsecs[0] << " rod2 " << acsecs[1]  
+	     << " rod3 " << acsecs[2] << G4endl;
+    }
+    
+    if (tot_abs_csec > small) {
+      G4double pw = -path * tot_abs_csec;
+      
+      if (pw < -huge_num) pw = -huge_num;
+      pw = 1.0 - std::exp(pw);
+      
+      if (verboseLevel > 2){
+	G4cout << " pw " << pw << G4endl;
+      }
+      
+      G4double apath = path;
+      
+      if (inuclRndm() < pw) 
+	apath = -1.0 / tot_abs_csec * std::log(1.0 - pw * inuclRndm());
+      
+      if (cparticle.young(young_cut, apath)) apath = path;  
+      
+      if(verboseLevel > 2){
+	G4cout << " apath " << apath << " path " << path << G4endl;
+      }
+      
+      if (apath < path) {	// choose the qdeutron
+	G4double sl = inuclRndm() * tot_abs_csec;
+	G4double as = 0.0;
+	
+	for (G4int i = 0; i < 3; i++) {
+	  as += acsecs[i];
+	  if (sl < as) { 
+	    if (verboseLevel > 2) G4cout << " deut type " << i << G4endl; 
+	    thePartners.push_back(partner(qdeutrons[i], apath));
+	    break;
+	  }
 	}
-
-	abs_sec = absorptionCrossSection(ekin, ptype); 
-	abs_sec *= nucleon_densities[1][zone] * nucleon_densities[1][zone] *
-	  rat1 * rat1 * vol; 
-
-      } else {
-	abs_sec = 0.0;
-      }; 
-
-      // abs_sec = 0.0;
-      tot_abs_csec += abs_sec;
-      acsecs.push_back(abs_sec);
-      qdeutrons.push_back(nnd);
-
-      if (verboseLevel > 2){
-	G4cout << " rod1 " << acsecs[0] << " rod2 " << acsecs[1]  
-	       << " rod3 " << acsecs[2] << G4endl;
-      }
-
-      if (tot_abs_csec > small) {
-     
-	G4double pw = -path * tot_abs_csec;
-
-	if (pw < -huge_num) pw = -huge_num;
-	pw = 1.0 - std::exp(pw);
-
-	if (verboseLevel > 2){
-	  G4cout << " pw " << pw << G4endl;
-	}
-
-	G4double apath = path;
-
-	if (inuclRndm() < pw) 
-	  apath = -1.0 / tot_abs_csec * std::log(1.0 - pw * inuclRndm());
-
-	if (cparticle.young(young_cut, apath)) apath = path;  
-
-	if(verboseLevel > 2){
-	  G4cout << " apath " << apath << " path " << path << G4endl;
-	}
-
-	if (apath < path) { // chose the qdeutron
-
-	  G4double sl = inuclRndm() * tot_abs_csec;
-	  G4double as = 0.0;
-
-	  for (G4int i = 0; i < 3; i++) {
-	    as += acsecs[i];
-
-	    if (sl < as) { 
-
-	      if (verboseLevel > 2){
-		G4cout << " deut type " << i << G4endl; 
-	      }
-
-	      thePartners.push_back(partner(qdeutrons[i], apath));
-
-	      break;
-	    };
-	  };
-	};    
-      };
-    };  
-
-    if(verboseLevel > 2){
-      G4cout << " after deutrons " << thePartners.size() << G4endl;
-    }
-  
-    if (thePartners.size() > 1) {		// Sort list by path length
-      std::sort(thePartners.begin(), thePartners.end(), sortPartners);
-    }
-
-    G4InuclElementaryParticle particle;		// Dummy for end of list
-    thePartners.push_back(partner(particle, path));
-  }
- 
-  return;
+      }    
+    }
+  }  
+  
+  if (verboseLevel > 2) {
+    G4cout << " after deuterons " << thePartners.size() << " partners"
+	   << G4endl;
+  }
+  
+  if (thePartners.size() > 1) {		// Sort list by path length
+    std::sort(thePartners.begin(), thePartners.end(), sortPartners);
+  }
+  
+  G4InuclElementaryParticle particle;		// Total path at end of list
+  thePartners.push_back(partner(particle, path));
 }
 
@@ -774,21 +796,30 @@
 G4NucleiModel::generateParticleFate(G4CascadParticle& cparticle,
                                     G4ElementaryParticleCollider* theElementaryParticleCollider) {
-  if (verboseLevel > 3)
+  if (verboseLevel > 1)
     G4cout << " >>> G4NucleiModel::generateParticleFate" << G4endl;
 
+  if (verboseLevel > 2) {
+    G4cout << " cparticle: ";
+    cparticle.print();
+  }
+
+  // Create four-vector checking
+#ifdef G4CASCADE_CHECK_ECONS
+  G4CascadeCheckBalance balance(0.005, 0.01, "G4NucleiModel");	// Second arg is in GeV
+  balance.setVerboseLevel(verboseLevel);
+#endif
+
   outgoing_cparticles.clear();		// Clear return buffer for this event
-
   generateInteractionPartners(cparticle);	// Fills "thePartners" data
 
   if (thePartners.empty()) { // smth. is wrong -> needs special treatment
-    G4cout << " generateParticleFate-> can not be here " << G4endl;
+    G4cerr << " generateParticleFate-> got empty interaction-partners list "
+	   << G4endl;
     return outgoing_cparticles;
   }
 
-  G4int npart = thePartners.size();
-
-  if (npart == 1) { // cparticle is on the next zone entry
-    // need to go here if particle outside nucleus ?
-    //
+  G4int npart = thePartners.size();	// Last item is a total-path placeholder
+
+  if (npart == 1) { 		// cparticle is on the next zone entry
     cparticle.propagateAlongThePath(thePartners[0].second);
     cparticle.incrementCurrentPath(thePartners[0].second);
@@ -796,54 +827,64 @@
     outgoing_cparticles.push_back(cparticle);
     
-    if (verboseLevel > 2){
+    if (verboseLevel > 2) {
       G4cout << " next zone " << G4endl;
       cparticle.print();
     }
-    
-  } else { // there are possible interactions
-    
+  } else {			// there are possible interactions
+    if (verboseLevel > 1)
+      G4cout << " processing " << npart-1 << " possible interactions" << G4endl;
+
     G4ThreeVector old_position = cparticle.getPosition();
-    
-    G4InuclElementaryParticle bullet = cparticle.getParticle();
-    
+    G4InuclElementaryParticle& bullet = cparticle.getParticle();
     G4bool no_interaction = true;
-    
     G4int zone = cparticle.getCurrentZone();
     
-    G4CollisionOutput output;
-
-    for (G4int i = 0; i < npart - 1; i++) {
+    for (G4int i=0; i<npart-1; i++) {	// Last item is a total-path placeholder
       if (i > 0) cparticle.updatePosition(old_position); 
       
-      G4InuclElementaryParticle target = thePartners[i].first; 
-      
-      if (verboseLevel > 2){
-	if (target.quasi_deutron()) 
-	  G4cout << " try absorption: target " << target.type() << " bullet " <<
-	    bullet.type() << G4endl;
+      G4InuclElementaryParticle& target = thePartners[i].first; 
+
+      if (verboseLevel > 3) {
+	if (target.quasi_deutron()) G4cout << " try absorption: ";
+	G4cout << " target " << target.type() << " bullet " << bullet.type()
+	       << G4endl;
       }
 
-      output.reset();
-      theElementaryParticleCollider->collide(&bullet, &target, output);
-      
-      if (verboseLevel > 2) output.printCollisionOutput();
+      EPCoutput.reset();
+      theElementaryParticleCollider->collide(&bullet, &target, EPCoutput);
+      
+      if (verboseLevel > 2) {
+	EPCoutput.printCollisionOutput();
+#ifdef G4CASCADE_CHECK_ECONS
+	balance.collide(&bullet, &target, EPCoutput);
+	balance.okay();		// Do checks, but ignore result
+#endif
+      }
 
       // Don't need to copy list, as "output" isn't changed again below
       const std::vector<G4InuclElementaryParticle>& outgoing_particles = 
-	output.getOutgoingParticles();
-      
-      if (passFermi(outgoing_particles, zone)) { // interaction
-	cparticle.propagateAlongThePath(thePartners[i].second);
-	G4ThreeVector new_position = cparticle.getPosition();
-	
-	for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 
-	  G4CascadParticle temp(outgoing_particles[ip], new_position, zone, 0.0, 0);
-	  outgoing_cparticles.push_back(temp);
-	}
-	
-	no_interaction = false;
-	current_nucl1 = 0;
-	current_nucl2 = 0;
-#ifdef CHC_CHECK
+	EPCoutput.getOutgoingParticles();
+      
+      if (!passFermi(outgoing_particles, zone)) continue; // Interaction fails
+
+      // Successful interaction, add results to output list
+      cparticle.propagateAlongThePath(thePartners[i].second);
+      G4ThreeVector new_position = cparticle.getPosition();
+
+      if (verboseLevel > 2)
+	G4cout << " adding " << outgoing_particles.size()
+	       << " output particles" << G4endl;
+      
+      for (G4int ip = 0; ip < G4int(outgoing_particles.size()); ip++) { 
+	G4CascadParticle temp(outgoing_particles[ip], new_position, zone, 0.0, 0);
+	outgoing_cparticles.push_back(temp);
+      }
+      
+      no_interaction = false;
+      current_nucl1 = 0;
+      current_nucl2 = 0;
+      
+#ifdef G4CASCADE_DEBUG_CHARGE
+      {
 	G4double out_charge = 0.0;
 	
@@ -851,50 +892,66 @@
 	  out_charge += outgoing_particles[ip].getCharge();
 	
-	G4cout << " multiplicity " << outgoing_particles.size() <<
-	  " bul type " << bullet.type() << " targ type " << target.type() << 
-	  G4endl << " initial charge " << bullet.getCharge() + target.getCharge() 
+	G4cout << " multiplicity " << outgoing_particles.size()
+	       << " bul type " << bullet.type()
+	       << " targ type " << target.type()
+	       << "\n initial charge "
+	       << bullet.getCharge() + target.getCharge() 
 	       << " out charge " << out_charge << G4endl;  
+      }
 #endif
-	
-	if (verboseLevel > 2){
-	  G4cout << " partner type " << target.type() << G4endl;
-	}
-	
-	if (target.nucleon()) {
-	  current_nucl1 = target.type();
-	  
-	} else {
-	  if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
-	  
-	  current_nucl1 = (target.type() - 100) / 10;
-	  current_nucl2 = target.type() - 100 - 10 * current_nucl1;
-	}   
-	
-	if (current_nucl1 == 1) {
-	  protonNumberCurrent -= 1.0;
-	  
-	} else {
-	  neutronNumberCurrent -= 1.0;
-	}; 
-	
-	if (current_nucl2 == 1) {
-	  protonNumberCurrent -= 1.0;
-	  
-	} else if(current_nucl2 == 2) {
-	  neutronNumberCurrent -= 1.0;
-	};
-	
-	break;
-      }; 
-    }  // loop over partners
-    
-    if (no_interaction) { // still no interactions
+      
+      if (verboseLevel > 2)
+	G4cout << " partner type " << target.type() << G4endl;
+      
+      if (target.nucleon()) {
+	current_nucl1 = target.type();
+      } else {
+	if (verboseLevel > 2) G4cout << " good absorption " << G4endl;
+	current_nucl1 = (target.type() - 100) / 10;
+	current_nucl2 = target.type() - 100 - 10 * current_nucl1;
+      }   
+      
+      if (current_nucl1 == 1) {
+	if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
+	protonNumberCurrent--;
+      } else {
+	if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
+	neutronNumberCurrent--;
+      } 
+      
+      if (current_nucl2 == 1) {
+	if (verboseLevel > 3) G4cout << " decrement proton count" << G4endl;
+	protonNumberCurrent--;
+      } else if (current_nucl2 == 2) {
+	if (verboseLevel > 3) G4cout << " decrement neutron count" << G4endl;
+	neutronNumberCurrent--;
+      }
+      
+      break;
+    }		// loop over partners
+    
+    if (no_interaction) { 		// still no interactions
+      if (verboseLevel > 1) G4cout << " no interaction " << G4endl;
+
+      // For conservation checking (below), get particle before updating
+      static G4InuclElementaryParticle prescatCP;	// Avoid memory churn
+      prescatCP = cparticle.getParticle();
+
+      // Last "partner" is just a total-path placeholder
       cparticle.updatePosition(old_position); 
-      cparticle.propagateAlongThePath(thePartners[npart - 1].second);
-      cparticle.incrementCurrentPath(thePartners[npart - 1].second);
+      cparticle.propagateAlongThePath(thePartners[npart-1].second);
+      cparticle.incrementCurrentPath(thePartners[npart-1].second);
       boundaryTransition(cparticle);
       outgoing_cparticles.push_back(cparticle);
-    };
-  }; 
+
+      // Check conservation for simple scattering (ignore target nucleus!)
+#ifdef G4CASCADE_CHECK_ECONS
+      if (verboseLevel > 2) {
+	balance.collide(&prescatCP, 0, outgoing_cparticles);
+	balance.okay();		// Report violations, but don't act on them
+      }
+#endif
+    }
+  }	// if (npart == 1) [else]
 
   return outgoing_cparticles;
@@ -903,91 +960,91 @@
 G4bool G4NucleiModel::passFermi(const std::vector<G4InuclElementaryParticle>& particles, 
 				G4int zone) {
+  if (verboseLevel > 1) {
+    G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
+  }
+
+  // Only check Fermi momenta for nucleons
+  for (G4int i = 0; i < G4int(particles.size()); i++) {
+    if (!particles[i].nucleon()) continue;
+
+    G4int type   = particles[i].type();
+    G4double mom = particles[i].getMomModule();
+    G4double pf  = fermi_momenta[type-1][zone];
+
+    if (verboseLevel > 2)
+      G4cout << " type " << type << " p " << mom << " pf " << pf << G4endl;
+    
+    if (mom < pf) {
+	if (verboseLevel > 2) G4cout << " rejected by Fermi" << G4endl;
+	return false;
+    }
+  }
+  return true; 
+}
+
+void G4NucleiModel::boundaryTransition(G4CascadParticle& cparticle) {
+  if (verboseLevel > 1) {
+    G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
+  }
+
+  G4int zone = cparticle.getCurrentZone();
+
+  if (cparticle.movingInsideNuclei() && zone == 0) {
+    G4cerr << " boundaryTransition-> in zone 0 " << G4endl;
+    return;
+  }
+
+  G4LorentzVector mom = cparticle.getMomentum();
+  G4ThreeVector pos = cparticle.getPosition();
+  
+  G4int type = cparticle.getParticle().type();
+  
+  G4double pr = pos.dot(mom.vect());
+  G4double r = pos.mag();
+  
+  pr /= r;
+  
+  G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
+  
+  G4double dv = getPotential(type,zone) - getPotential(type, next_zone);
+  //    G4cout << "Potentials for type " << type << " = " 
+  //           << getPotential(type,zone) << " , "
+  //	   << getPotential(type,next_zone) << G4endl;
+  
+  G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr;
+  
+  G4double p1r;
+  
   if (verboseLevel > 3) {
-    G4cout << " >>> G4NucleiModel::passFermi" << G4endl;
-  }
-
-  for (G4int i = 0; i < G4int(particles.size()); i++) {
-
-    if (particles[i].nucleon()) {
-
-      if (verboseLevel > 2){
-	G4cout << " type " << particles[i].type() << " p " << particles[i].getMomModule()
-	       << " pf " << fermi_momenta[particles[i].type() - 1][zone] << G4endl;
-      }
-
-      if (particles[i].getMomModule() < fermi_momenta[particles[i].type() - 1][zone]) {
-
-	if (verboseLevel > 2) {
-	  G4cout << " rejected by fermi: type " << particles[i].type() << 
-	    " p " << particles[i].getMomModule() << G4endl;
-	}
-
-	return false;
-      };
-    };
-  };
-  return true; 
-}
-
-void G4NucleiModel::boundaryTransition(G4CascadParticle& cparticle) {
-
+    G4cout << " type " << type << " zone " << zone << " next " << next_zone
+	   << " qv " << qv << " dv " << dv << G4endl;
+  }
+
+  if (qv <= 0.0) { 	// reflection
+    if (verboseLevel > 3) G4cout << " reflects off boundary" << G4endl;
+    p1r = -pr;
+    cparticle.incrementReflectionCounter();
+  } else {		// transition
+    if (verboseLevel > 3) G4cout << " passes thru boundary" << G4endl;
+    p1r = std::sqrt(qv);
+    if(pr < 0.0) p1r = -p1r;
+    cparticle.updateZone(next_zone);
+    cparticle.resetReflection();
+  }
+  
+  G4double prr = (p1r - pr) / r;  
+  
   if (verboseLevel > 3) {
-    G4cout << " >>> G4NucleiModel::boundaryTransition" << G4endl;
-  }
-
-  G4int zone = cparticle.getCurrentZone();
-
-  if (cparticle.movingInsideNuclei() && zone == 0) {
-    G4cout << " boundaryTransition-> in zone 0 " << G4endl;
-
-  } else {
-    G4LorentzVector mom = cparticle.getMomentum();
-    G4ThreeVector pos = cparticle.getPosition();
-
-    G4int type = cparticle.getParticle().type();
-
-    G4double pr = pos.dot(mom.vect());
-    G4double r = pos.mag();
-
-    pr /= r;
-
-    G4int next_zone = cparticle.movingInsideNuclei() ? zone - 1 : zone + 1;
-
-    G4double dv = getPotential(type,zone) - getPotential(type, next_zone);
-    //    G4cout << "Potentials for type " << type << " = " 
-    //           << getPotential(type,zone) << " , "
-    //	   << getPotential(type,next_zone) << G4endl;
-
-    G4double qv = dv * dv - 2.0 * dv * mom.e() + pr * pr;
-
-    G4double p1r;
-
-    if (verboseLevel > 2){
-      G4cout << " type " << type << " zone " << zone
-             << " next " << next_zone
-             << " qv " << qv << " dv " << dv << G4endl;
-    }
-
-    if(qv <= 0.0) { // reflection 
-      p1r = -pr;
-      cparticle.incrementReflectionCounter();
-
-    } else { // transition
-      p1r = std::sqrt(qv);
-      if(pr < 0.0) p1r = -p1r;
-      cparticle.updateZone(next_zone);
-      cparticle.resetReflection();
-    };
- 
-    G4double prr = (p1r - pr) / r;  
-
-    mom.setVect(mom.vect() + pos*prr);
-    cparticle.updateParticleMomentum(mom);
-  }; 
+    G4cout << " prr " << prr << " delta px " << prr*pos.x() << " py "
+	   << prr*pos.y()  << " pz " << prr*pos.z() << " mag "
+	   << std::fabs(prr*r) << G4endl;
+  }
+
+  mom.setVect(mom.vect() + pos*prr);
+  cparticle.updateParticleMomentum(mom);
 }
 
 G4bool G4NucleiModel::worthToPropagate(const G4CascadParticle& cparticle) const {
-
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::worthToPropagate" << G4endl;
   }
@@ -1010,69 +1067,48 @@
 	     << worth << G4endl;
     }
-  };
+  }
 
   return worth;
 }
 
+
 G4double G4NucleiModel::getRatio(G4int ip) const {
-
-  if (verboseLevel > 3) {
+  if (verboseLevel > 1) {
     G4cout << " >>> G4NucleiModel::getRatio" << G4endl;
   }
 
-  G4double rat;
-  //  G4double ratm;
-
-  // Calculate number of protons and neutrons in local region
-  //  G4double Athird = G4cbrt(A);
-  //  G4double Nneut = Athird*(A-Z)/A;
-  //  G4double Nprot = Athird*Z/A;
-
-  // Reduce number of 
   if (ip == 1) {
+    if (verboseLevel > 2) {
+      G4cout << " current " << protonNumberCurrent << " inp " << protonNumber
+	     << G4endl;
+    }
+    return G4double(protonNumberCurrent)/G4double(protonNumber);
+  }
+
+  if (ip == 2) {
     if (verboseLevel > 2){
-      G4cout << " current " << protonNumberCurrent << " inp " << protonNumber << G4endl;
-    }
-
-    rat = protonNumberCurrent/protonNumber;
-
-    // Calculate ratio modified for local region
-    //    G4double deltaP = protonNumber - protonNumberCurrent;
-    //    G4cout << " deltaP = " << deltaP << G4endl;
-    //    ratm = std::max(0.0, (Nprot - deltaP)/Nprot);
-
-  } else {
-    if (verboseLevel > 2){
-      G4cout << " current " << neutronNumberCurrent << " inp " << neutronNumber << G4endl;
-    }
-
-    rat = neutronNumberCurrent/neutronNumber;
-
-    // Calculate ratio modified for local region
-    //    G4double deltaN = neutronNumber - neutronNumberCurrent;
-    //   G4cout << " deltaN = " << deltaN << G4endl;
-    //    ratm = std::max(0.0, (Nneut - deltaN)/Nneut);
-  }
-
-  //  G4cout << " get ratio: ratm =  " << ratm << G4endl;
-  return rat;
-  //  return ratm;
-}
-
-G4CascadParticle G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle) {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle)" << G4endl;
+      G4cout << " current " << neutronNumberCurrent << " inp " << neutronNumber
+	     << G4endl;
+    }
+    return G4double(neutronNumberCurrent)/G4double(neutronNumber);
+  }
+
+  return 0.;
+}
+
+G4CascadParticle 
+G4NucleiModel::initializeCascad(G4InuclElementaryParticle* particle) {
+  if (verboseLevel > 1) {
+    G4cout << " >>> G4NucleiModel::initializeCascad(particle)" << G4endl;
   }
 
   const G4double large = 1000.0;
 
-  G4double s1 = std::sqrt(inuclRndm()); 
-  G4double phi = randomPHI();
-  G4double rz = nuclei_radius * s1;
-
-  G4ThreeVector pos(rz*std::cos(phi), rz*std::sin(phi),
-		    -nuclei_radius*std::sqrt(1.0 - s1*s1));
- 
+  // FIXME:  Previous version generated random sin(theta), then used -cos(theta)
+  //         Using generateWithRandomAngles changes result!
+  // G4ThreeVector pos = generateWithRandomAngles(nuclei_radius).vect();
+  G4double costh = std::sqrt(1.0 - inuclRndm());
+  G4ThreeVector pos = generateWithFixedTheta(-costh, nuclei_radius);
+
   G4CascadParticle cpart(*particle, pos, number_of_zones, large, 0);
 
@@ -1085,7 +1121,7 @@
 				     G4InuclNuclei* target,
 				     modelLists& output) {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4NucleiModel::initializeCascad(G4InuclNuclei* bullet, G4InuclNuclei* target)" << G4endl;
+  if (verboseLevel) {
+    G4cout << " >>> G4NucleiModel::initializeCascad(bullet,target,output)"
+	   << G4endl;
   }
 
@@ -1110,8 +1146,8 @@
 
   // first decide whether it will be cascad or compound final nuclei
-  G4double ab = bullet->getA();
-  G4double zb = bullet->getZ();
-  G4double at = target->getA();
-  G4double zt = target->getZ();
+  G4int ab = bullet->getA();
+  G4int zb = bullet->getZ();
+  G4int at = target->getA();
+  G4int zt = target->getZ();
 
   G4double massb = bullet->getMass();	// For creating LorentzVectors below
@@ -1119,6 +1155,6 @@
   if (ab < max_a_for_cascad) {
 
-    G4double benb = 0.001 * G4NucleiProperties::GetBindingEnergy(G4lrint(ab), G4lrint(zb)) / ab;
-    G4double bent = 0.001 * G4NucleiProperties::GetBindingEnergy(G4lrint(at), G4lrint(zt)) / at;
+    G4double benb = 0.001 * bindingEnergy(ab,zb) / G4double(ab);
+    G4double bent = 0.001 * bindingEnergy(at,zt) / G4double(at);
     G4double ben = benb < bent ? bent : benb;
 
@@ -1128,22 +1164,15 @@
       while (casparticles.size() == 0 && itryg < itry_max) {      
 	itryg++;
-
-	if(itryg > 0) particles.clear();
+	particles.clear();
       
 	//    nucleons coordinates and momenta in nuclei rest frame
-	std::vector<G4ThreeVector> coordinates;
-	std::vector<G4LorentzVector> momentums;
+	coordinates.clear();
+	momentums.clear();
      
-	if (ab < 3.0) { // deutron, simplest case
+	if (ab < 3) { // deuteron, simplest case
 	  G4double r = 2.214 - 3.4208 * std::log(1.0 - 0.981 * inuclRndm());
-	  G4double s = 2.0 * inuclRndm() - 1.0;
-	  G4double r1 = r * std::sqrt(1.0 - s * s);
-	  G4double phi = randomPHI();
-
-	  G4ThreeVector coord1(r1*std::cos(phi), r1*std::sin(phi), r*s);   
+	  G4ThreeVector coord1 = generateWithRandomAngles(r).vect();
 	  coordinates.push_back(coord1);
-
-	  coord1 *= -1.;
-	  coordinates.push_back(coord1);
+	  coordinates.push_back(-coord1);
 
 	  G4double p = 0.0;
@@ -1157,5 +1186,5 @@
 	    if (p * p / (p * p + 2079.36) / (p * p + 2079.36) > 1.2023e-4 * inuclRndm() &&
 	       p * r > 312.0) bad = false;
-	  };
+	  }
 
 	  if (itry == itry_max)
@@ -1176,6 +1205,4 @@
 	  momentums.push_back(-mom);
 	} else {
-	  G4int ia = int(ab + 0.5);
-
 	  G4ThreeVector coord1;
 
@@ -1184,5 +1211,5 @@
 	  G4int itry = 0;
         
-	  if (ab < 4.0) { // a == 3
+	  if (ab == 3) {
 	    while (badco && itry < itry_max) {
 	      if (itry > 0) coordinates.clear();
@@ -1210,6 +1237,6 @@
 		    }
 		    break;
-		  };
-		};
+		  }
+		}
 
 		if (itry1 == itry_max) { // bad case
@@ -1217,6 +1244,6 @@
 		  coordinates.push_back(coord1);
 		  break;
-		};
-	      };
+		}
+	      }
 
 	      coord1 = -coordinates[0] - coordinates[1]; 
@@ -1241,13 +1268,13 @@
 
 		    break; 
-		  };      
-		};
+		  }      
+		}
 
 		if (large_dist) break;
-	      }; 
+	      } 
 
 	      if(!large_dist) badco = false;
 
-	    };
+	    }
 
 	  } else { // a >= 4
@@ -1264,5 +1291,5 @@
 	      G4int i(0);
 	    
-	      for (i = 0; i < ia-1; i++) {
+	      for (i = 0; i < ab-1; i++) {
 		G4int itry1 = 0;
 		G4double s, u; 
@@ -1283,6 +1310,6 @@
 
 		    break;
-		  };
-		};
+		  }
+		}
 
 		if (itry1 == itry_max) { // bad case
@@ -1290,9 +1317,9 @@
 		  coordinates.push_back(coord1);
 		  break;
-		};
-	      };
+		}
+	      }
 
 	      coord1 *= 0.0;	// Cheap way to reset
-	      for(G4int j = 0; j < ia -1; j++) coord1 -= coordinates[j];
+	      for(G4int j = 0; j < ab -1; j++) coord1 -= coordinates[j];
 
 	      coordinates.push_back(coord1);   
@@ -1304,6 +1331,6 @@
 	      G4bool large_dist = false;
 
-	      for (i = 0; i < ia-1; i++) {
-		for (G4int j = i+1; j < ia; j++) {
+	      for (i = 0; i < ab-1; i++) {
+		for (G4int j = i+1; j < ab; j++) {
 	     
 		  G4double r2 = (coordinates[i]-coordinates[j]).mag2();
@@ -1317,13 +1344,13 @@
 
 		    break; 
-		  };      
-		};
+		  }      
+		}
 
 		if (large_dist) break;
-	      }; 
+	      } 
 
 	      if (!large_dist) badco = false;
-	    };
-	  }; 
+	    }
+	  } 
 
 	  if(badco) {
@@ -1339,7 +1366,6 @@
 	    G4LorentzVector mom;
 	    //G4bool badp = True;
-	    G4int i(0);
-
-	    for (i = 0; i < ia - 1; i++) {
+
+	    for (G4int i = 0; i < ab - 1; i++) {
 	      G4int itry = 0;
 
@@ -1355,6 +1381,6 @@
 
 		  break;
-		};
-	      };
+		}
+	      }
 
 	      if(itry == itry_max) {
@@ -1365,14 +1391,16 @@
 		particles.clear();
 		return;
-	      }; 
-
-	    };
+	      } 
+
+	    }
 	    // last momentum
 
 	    mom *= 0.;		// Cheap way to reset
-	    for(G4int j=0; j< ia-1; j++) mom -= momentums[j]; 
+	    mom.setE(bullet->getEnergy()+target->getEnergy());
+
+	    for(G4int j=0; j< ab-1; j++) mom -= momentums[j]; 
 
 	    momentums.push_back(mom);
-	  }; 
+	  } 
 	}
  
@@ -1385,5 +1413,5 @@
 
 	  if(rp > rb) rb = rp;
-	};
+	}
 
 	// nuclei i.p. as a whole
@@ -1396,20 +1424,16 @@
 	for (i = 0; i < G4int(coordinates.size()); i++) {
 	  coordinates[i] += global_pos;
-	};  
+	}  
 
 	// all nucleons at rest
-	std::vector<G4InuclElementaryParticle> raw_particles;
-	G4int ia = G4int(ab + 0.5);
-	G4int iz = G4int(zb + 0.5);
-
-	for (G4int ipa = 0; ipa < ia; ipa++) {
-	  G4int knd = ipa < iz ? 1 : 2;
+	raw_particles.clear();
+
+	for (G4int ipa = 0; ipa < ab; ipa++) {
+	  G4int knd = ipa < zb ? 1 : 2;
 	  raw_particles.push_back(G4InuclElementaryParticle(momentums[ipa], knd));
-	}; 
+	} 
       
 	G4InuclElementaryParticle dummy(small_ekin, 1);
-	G4LorentzConvertor toTheBulletRestFrame;
-	toTheBulletRestFrame.setBullet(dummy.getMomentum(), dummy.getMass());
-	toTheBulletRestFrame.setTarget(bullet->getMomentum(),bullet->getMass());
+	G4LorentzConvertor toTheBulletRestFrame(&dummy, bullet);
 	toTheBulletRestFrame.toTheTargetRestFrame();
 
@@ -1418,5 +1442,5 @@
 	for (ipart = raw_particles.begin(); ipart != raw_particles.end(); ipart++) {
 	  ipart->setMomentum(toTheBulletRestFrame.backToTheLab(ipart->getMomentum())); 
-	};
+	}
 
 	// fill cascad particles and outgoing particles
@@ -1437,10 +1461,10 @@
 	      if(t1 > 0.0) {
 		if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
-	      };
+	      }
 
 	      if(tr < 0.0 && t2 > 0.0) {
 
 		if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
-	      };
+	      }
 
 	    } else {
@@ -1448,15 +1472,15 @@
 
 		if(coordinates[ip].z() + mom.z() * t2 / pmod <= 0.0) tr = t2;
-	      };
+	      }
 
 	      if(tr < 0.0 && t1 > 0.0) {
 		if(coordinates[ip].z() + mom.z() * t1 / pmod <= 0.0) tr = t1;
-	      };
-	    }; 
-
-	  };
+	      }
+	    } 
+
+	  }
 
 	  if(tr >= 0.0) { // cascad particle
-	    coordinates[ip] += mom*tr / pmod;
+	    coordinates[ip] += mom.vect()*tr / pmod;
 	    casparticles.push_back(G4CascadParticle(raw_particles[ip], 
 						    coordinates[ip], 
@@ -1465,7 +1489,7 @@
 	  } else {
 	    particles.push_back(raw_particles[ip]); 
-	  }; 
-	};
-      };    
+	  } 
+	}
+      }    
 
       if(casparticles.size() == 0) {
@@ -1474,7 +1498,7 @@
 	G4cout << " can not generate proper distribution for " << itry_max
 	       << " steps " << G4endl;
-      };    
-    };
-  };
+      }    
+    }
+  }
 
   if(verboseLevel > 2){
@@ -1509,5 +1533,5 @@
   } else if (ke < 1.0) {
     csec = 3.6735 * (1.0-ke)*(1.0-ke);     
-  };
+  }
 
   if (csec < 0.0) csec = 0.0;
