Index: trunk/source/processes/hadronic/models/cascade/cascade/src/G4IntraNucleiCascader.cc
===================================================================
--- trunk/source/processes/hadronic/models/cascade/cascade/src/G4IntraNucleiCascader.cc	(revision 1007)
+++ trunk/source/processes/hadronic/models/cascade/cascade/src/G4IntraNucleiCascader.cc	(revision 1315)
@@ -23,34 +23,53 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
+// $Id: G4IntraNucleiCascader.cc,v 1.35 2010/05/21 17:56:34 mkelsey Exp $
+// Geant4 tag: $Name: geant4-09-04-beta-cand-01 $
 //
-
+// 20100114  M. Kelsey -- Remove G4CascadeMomentum, use G4LorentzVector directly
+// 20100307  M. Kelsey -- Bug fix: momentum_out[0] should be momentum_out.e()
+// 20100309  M. Kelsey -- Eliminate some unnecessary std::pow()
+// 20100407  M. Kelsey -- Pass "all_particles" as argument to initializeCascad,
+//		following recent change to G4NucleiModel.
+// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
+// 20100517  M. Kelsey -- Inherit from common base class, make other colliders
+//		simple data members
+
+#include "G4IntraNucleiCascader.hh"
 #define RUN
 
-#include "G4IntraNucleiCascader.hh"
+#include "G4CascadParticle.hh"
+#include "G4CollisionOutput.hh"
+#include "G4ElementaryParticleCollider.hh"
+#include "G4HadTmpUtil.hh"
 #include "G4InuclElementaryParticle.hh"
 #include "G4InuclNuclei.hh"
-#include "G4LorentzConvertor.hh"
+#include "G4InuclSpecialFunctions.hh"
+#include "G4NucleiModel.hh"
+#include "G4NucleiProperties.hh"
 #include "G4ParticleLargerEkin.hh"
-#include "G4NucleiModel.hh"
-#include "G4CascadParticle.hh"
 #include "Randomize.hh"
 #include <algorithm>
 
+using namespace G4InuclSpecialFunctions;
+
+
 typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
 
 G4IntraNucleiCascader::G4IntraNucleiCascader()
-  : verboseLevel(1) {
+  : G4VCascadeCollider("G4IntraNucleiCascader"),
+    theElementaryParticleCollider(new G4ElementaryParticleCollider) {}
+
+G4IntraNucleiCascader::~G4IntraNucleiCascader() {
+  delete theElementaryParticleCollider;
+}
+
+
+void G4IntraNucleiCascader::collide(G4InuclParticle* bullet,
+				    G4InuclParticle* target,
+				    G4CollisionOutput& output) {
 
   if (verboseLevel > 3) {
-    G4cout << " >>> G4IntraNucleiCascader::G4IntraNucleiCascader" << G4endl;
-  }
-
-}
-
-G4CollisionOutput G4IntraNucleiCascader::collide(G4InuclParticle* bullet,
-						 G4InuclParticle* target) {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4IntraNucleiCascader::collide" << G4endl;
+    G4cout << " >>> G4IntraNucleiCascader::collide inter_case " << inter_case 
+	   << G4endl;
   }
 
@@ -63,6 +82,4 @@
     target->printParticle();
   }
-
-  G4CollisionOutput output;
 
 #ifdef RUN
@@ -73,9 +90,9 @@
   G4NucleiModel model(tnuclei);
   G4double coulombBarrier = 0.00126*tnuclei->getZ()/
-                                      (1.+std::pow(tnuclei->getA(),0.333));
-
-  G4CascadeMomentum momentum_in = bullet->getMomentum();
-
-  momentum_in[0] += tnuclei->getMass();
+                                      (1.+G4cbrt(tnuclei->getA()));
+
+  G4LorentzVector momentum_in = bullet->getMomentum();
+
+  momentum_in.setE(momentum_in.e()+tnuclei->getMass());
 
   G4double ekin_in; 
@@ -83,6 +100,7 @@
   if (verboseLevel > 3) {
     model.printModel();
-    G4cout << " intitial momentum  E " << momentum_in[0] << " Px " << momentum_in[1] 
-	   << " Py " << momentum_in[2] << " Pz " << momentum_in[3] << G4endl;
+    G4cout << " intitial momentum  E " << momentum_in.e() << " Px "
+	   << momentum_in.x() << " Py " << momentum_in.y() << " Pz "
+	   << momentum_in.z() << G4endl;
   }
 
@@ -117,6 +135,6 @@
       zfin += zb;
 
-      std::pair<std::vector<G4CascadParticle>, std::vector<G4InuclElementaryParticle> > 
-	all_particles = model.initializeCascad(bnuclei, tnuclei);
+      G4NucleiModel::modelLists all_particles;    // Buffer to receive lists
+      model.initializeCascad(bnuclei, tnuclei, all_particles);
 
       cascad_particles = all_particles.first;
@@ -139,5 +157,4 @@
 
 	for (i = 0; i < ihn; i++) theExitonConfiguration.incrementHoles(2);
-
 	for (i = 0; i < ihz; i++) theExitonConfiguration.incrementHoles(1);
       };
@@ -208,4 +225,5 @@
                                             (1./KE - 1./coulombBarrier)*
                                          std::sqrt(mass*(coulombBarrier-KE)) );
+
             if (G4UniformRand() < CBP) {
 	      output_particles.push_back(currentParticle);
@@ -241,11 +259,9 @@
     }
 
-    G4CascadeMomentum momentum_out;
+    G4LorentzVector momentum_out;
     particleIterator ipart;
 
     for (ipart = output_particles.begin(); ipart != output_particles.end(); ipart++) {
-      const G4CascadeMomentum& mom = ipart->getMomentum();
-
-      for (G4int j = 0; j < 4; j++) momentum_out[j] += mom[j];
+      momentum_out += ipart->getMomentum();
 
       zfin -= ipart->getCharge();
@@ -260,18 +276,16 @@
       G4InuclNuclei outgoing_nuclei(afin, zfin);
       G4double mass = outgoing_nuclei.getMass();
-      momentum_out[0] += mass;        
-
-      for (int j = 0; j < 4; j++) momentum_out[j] = momentum_in[j] - momentum_out[j];
+      momentum_out.setE(momentum_out.e()+mass);
+
+      momentum_out = momentum_in - momentum_out;
 
       if (verboseLevel > 3) {
-	G4cout << "  Eex + Ekin " << momentum_out[0]  <<  G4endl;
+	G4cout << "  Eex + Ekin " << momentum_out.e()  <<  G4endl;
       }
 
-      if (momentum_out[0] > 0.0) { // Eex + Ekin > 0.0
-	G4double pnuc = momentum_out[1] * momentum_out[1] + 
-	  momentum_out[2] * momentum_out[2] +
-	  momentum_out[3] * momentum_out[3]; 
+      if (momentum_out.e() > 0.0) { // Eex + Ekin > 0.0
+	G4double pnuc = momentum_out.vect().mag2(); 
 	G4double ekin = std::sqrt(mass * mass + pnuc) - mass;
-	G4double Eex = 1000.0 * (momentum_out[0] - ekin);
+	G4double Eex = 1000.0 * (momentum_out.e() - ekin);
 
 	if (verboseLevel > 3) {
@@ -283,10 +297,9 @@
 	  output.addOutgoingParticles(output_particles);
 	  outgoing_nuclei.setMomentum(momentum_out);
-	  outgoing_nuclei.setEnergy();
 	  outgoing_nuclei.setExitationEnergy(Eex);
 	  outgoing_nuclei.setExitonConfiguration(theExitonConfiguration);	                           	  
           output.addTargetFragment(outgoing_nuclei);
 
-	  return output;
+	  return;
 	};
       };
@@ -296,5 +309,5 @@
       if (afin == 1.0) { // recoiling nucleon
 
-	for (int j = 0; j < 4; j++) momentum_out[j] = momentum_in[j] - momentum_out[j];
+	momentum_out = momentum_in - momentum_out;
 
 	G4InuclElementaryParticle  last_particle;
@@ -314,5 +327,5 @@
       output.addOutgoingParticles(output_particles);
 
-      return output;
+      return;
     }; 
   };
@@ -322,9 +335,8 @@
   // special branch to avoid the cascad generation but to get the input for evaporation etc
 
-  G4CascadeMomentum momentum_out;
+  G4LorentzVector momentum_out;
   G4InuclNuclei outgoing_nuclei(169, 69);
 
   outgoing_nuclei.setMomentum(momentum_out);
-  outgoing_nuclei.setEnergy();
   outgoing_nuclei.setExitationEnergy(150.0);
 
@@ -334,5 +346,5 @@
   output.addTargetFragment(outgoing_nuclei);
 
-  return output;
+  return;
 
   /*
@@ -342,5 +354,5 @@
     output.addOutgoingParticle(*bparticle);
     output.addTargetFragment(*tnuclei);
-    return output;
+    return;
   */
 
@@ -354,5 +366,5 @@
   output.trivialise(bullet, target);
 
-  return output;
+  return;
 }
 
@@ -374,5 +386,6 @@
   if (eexs > eexs_cut) {
     G4double eexs_max0z = 1000.0 * ein / ediv_cut;
-    G4double dm = bindingEnergy(a, z);
+    //    G4double dm = bindingEnergy(a, z);
+    G4double dm = G4NucleiProperties::GetBindingEnergy(G4lrint(a), G4lrint(z));
     G4double eexs_max = eexs_max0z > reason_cut*dm ? eexs_max0z : reason_cut * dm;
 
