Index: trunk/source/processes/hadronic/models/cascade/cascade/src/G4InuclCollider.cc
===================================================================
--- trunk/source/processes/hadronic/models/cascade/cascade/src/G4InuclCollider.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/cascade/cascade/src/G4InuclCollider.cc	(revision 1315)
@@ -23,26 +23,45 @@
 // * acceptance of all terms of the Geant4 Software license.          *
 // ********************************************************************
+// $Id: G4InuclCollider.cc,v 1.30 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
+// 20100309  M. Kelsey -- Eliminate some unnecessary std::pow()
+// 20100413  M. Kelsey -- Pass G4CollisionOutput by ref to ::collide()
+// 20100418  M. Kelsey -- Move lab-frame transformation code to G4CollisonOutput
+// 20100429  M. Kelsey -- Change "photon()" to "isPhoton()"
+// 20100517  M. Kelsey -- Inherit from common base class, make other colliders
+//		simple data members, consolidate code
+
 #include "G4InuclCollider.hh"
+#include "G4BigBanger.hh"
+#include "G4CollisionOutput.hh"
+#include "G4ElementaryParticleCollider.hh"
+#include "G4EquilibriumEvaporator.hh"
+#include "G4IntraNucleiCascader.hh"
 #include "G4InuclElementaryParticle.hh"
 #include "G4LorentzConvertor.hh"
-#include "G4ParticleLargerEkin.hh"
-#include <algorithm>
-
-typedef std::vector<G4InuclElementaryParticle>::iterator particleIterator;
-typedef std::vector<G4InuclNuclei>::iterator nucleiIterator;
-	 
+#include "G4NonEquilibriumEvaporator.hh"
+
+
 G4InuclCollider::G4InuclCollider()
-  : verboseLevel(0) {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4InuclCollider::G4InuclCollider" << G4endl;
-  }
+  : G4VCascadeCollider("G4InuclCollider"),
+    theElementaryParticleCollider(new G4ElementaryParticleCollider),
+    theIntraNucleiCascader(new G4IntraNucleiCascader),
+    theNonEquilibriumEvaporator(new G4NonEquilibriumEvaporator),
+    theEquilibriumEvaporator(new G4EquilibriumEvaporator),
+    theBigBanger(new G4BigBanger) {}
+
+G4InuclCollider::~G4InuclCollider() {
+  delete theElementaryParticleCollider;
+  delete theIntraNucleiCascader;
+  delete theNonEquilibriumEvaporator;
+  delete theEquilibriumEvaporator;
+  delete theBigBanger;
 }
 
-G4CollisionOutput G4InuclCollider::collide(G4InuclParticle* bullet,
-					   G4InuclParticle* target) {
-
-  verboseLevel = 0;
+
+void G4InuclCollider::collide(G4InuclParticle* bullet, G4InuclParticle* target,
+			      G4CollisionOutput& globalOutput) {
   if (verboseLevel > 3) {
     G4cout << " >>> G4InuclCollider::collide" << G4endl;
@@ -51,29 +70,20 @@
   const G4int itry_max = 1000;
   		     
-  G4CollisionOutput globalOutput;
-  G4InuclElementaryParticle* particle1 =
-    dynamic_cast<G4InuclElementaryParticle*>(bullet);
-  G4InuclElementaryParticle* particle2 =
-    dynamic_cast<G4InuclElementaryParticle*>(target);
-  
-  if (particle1 && particle2) { // particle + particle (NOTE: also the h + H(1,1) treated here) 
+  if (useEPCollider(bullet,target)) {
     if (verboseLevel > 2) {
-      particle1->printParticle();
-      particle2->printParticle();
+      bullet->printParticle();
+      target->printParticle();
     }
 
-    globalOutput = theElementaryParticleCollider->collide(bullet, target);
-
+    theElementaryParticleCollider->collide(bullet, target, globalOutput);
   } else { // needs to call all machinery    	
     G4LorentzConvertor convertToTargetRestFrame;
-    G4InteractionCase interCase = bulletTargetSetter(bullet, target);
-    G4int intcase = interCase.getInterCase();
-     
-    if (intcase > 0) { // ok
+
+    interCase.set(bullet,target);
+    if (interCase.valid()) { // ok
       G4InuclNuclei* ntarget =
 	dynamic_cast<G4InuclNuclei*>(interCase.getTarget());
 
-      convertToTargetRestFrame.setTarget(ntarget->getMomentum(),
-					 ntarget->getMass());
+      convertToTargetRestFrame.setTarget(ntarget);
       G4int btype = 0;
       G4double ab = 0.0;
@@ -82,17 +92,15 @@
       G4double zt = ntarget->getZ();
        
-      if (intcase == 1) { // particle with nuclei
+      if (interCase.hadNucleus()) { // particle with nuclei
 	G4InuclElementaryParticle* pbullet = 
 	  dynamic_cast<G4InuclElementaryParticle*>(interCase.getBullet());
          
-	if (pbullet->photon()) {
-	  G4cout << " InuclCollider -> can not collide with photon " << G4endl;
+	if (pbullet->isPhoton()) {
+	  G4cerr << " InuclCollider -> can not collide with photon " << G4endl;
 
 	  globalOutput.trivialise(bullet, target);
-
-	  return globalOutput;
+	  return;
 	} else {
-	  convertToTargetRestFrame.setBullet(pbullet->getMomentum(),
-					     pbullet->getMass());   
+	  convertToTargetRestFrame.setBullet(pbullet);   
 	  btype = pbullet->type();
 	}; 
@@ -102,6 +110,5 @@
 	  dynamic_cast<G4InuclNuclei*>(interCase.getBullet());
 
-	convertToTargetRestFrame.setBullet(nbullet->getMomentum(),
-					   nbullet->getMass());   
+	convertToTargetRestFrame.setBullet(nbullet);   
 	ab = nbullet->getA();
 	zb = nbullet->getZ();
@@ -121,44 +128,41 @@
 	}
 
-	G4CascadeMomentum bmom;
-
-	bmom[3] = convertToTargetRestFrame.getTRSMomentum();
-
-	G4InuclNuclei ntarget(at, zt);
-	G4CascadeMomentum tmom;
-
-	ntarget.setMomentum(tmom);
-	ntarget.setEnergy();
-	theIntraNucleiCascader->setInteractionCase(intcase);
+	G4LorentzVector bmom;
+	bmom.setZ(convertToTargetRestFrame.getTRSMomentum());
+
+	G4InuclNuclei ntarget(at, zt);		// Default is at rest
+
+	theIntraNucleiCascader->setInteractionCase(interCase.code());
 	 
 	G4bool bad = true;
 	G4int itry = 0;
 	 
+	G4CollisionOutput TRFoutput;
+	G4CollisionOutput output;
 	while (bad && itry < itry_max) {
-	  G4CollisionOutput TRFoutput;
-	  G4CollisionOutput output;
-
 	  itry++;
-	  if (intcase == 1) {
+
+	  output.reset();	// Clear buffers for this attempt
+	  TRFoutput.reset();
+
+	  if (interCase.hadNucleus()) {
 	    G4InuclElementaryParticle pbullet(bmom, btype);
 
-	    output = theIntraNucleiCascader->collide(&pbullet, &ntarget);
+	    theIntraNucleiCascader->collide(&pbullet, &ntarget, output);
 	  } else {
-	    G4InuclNuclei nbullet(ab, zb);
-	    nbullet.setMomentum(bmom);
-	    nbullet.setEnergy();
-	    output = theIntraNucleiCascader->collide(&nbullet, &ntarget);
+	    G4InuclNuclei nbullet(bmom, ab, zb);
+	    theIntraNucleiCascader->collide(&nbullet, &ntarget, output);
 	  };   
 
 	  if (verboseLevel > 3) {
 	    G4cout << " After Cascade " << G4endl;
-
 	    output.printCollisionOutput();
 	  }
 	  
 	  // the rest, if any
+	  // FIXME:  The code below still does too much copying!
 	  TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
 
-	  if (output.numberOfNucleiFragments() == 1) { // there is smth. after	    
+	  if (output.numberOfNucleiFragments() == 1) { // there is smth. after
 	    G4InuclNuclei cascad_rec_nuclei = output.getNucleiFragments()[0];
 	    if (explosion(&cascad_rec_nuclei)) {
@@ -167,8 +171,8 @@
 	      };
 
-	      output = theBigBanger->collide(0,&cascad_rec_nuclei);
-	      TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
+	      theBigBanger->collide(0,&cascad_rec_nuclei, TRFoutput);
 	    } else {
-	      output = theNonEquilibriumEvaporator->collide(0, &cascad_rec_nuclei);
+	      output.reset();
+	      theNonEquilibriumEvaporator->collide(0, &cascad_rec_nuclei, output);
 
 	      if (verboseLevel > 3) {
@@ -179,62 +183,27 @@
 	      TRFoutput.addOutgoingParticles(output.getOutgoingParticles());
 	      G4InuclNuclei exiton_rec_nuclei = output.getNucleiFragments()[0];
-	      output = theEquilibriumEvaporator->collide(0, &exiton_rec_nuclei);
+
+	      output.reset();
+	      theEquilibriumEvaporator->collide(0, &exiton_rec_nuclei, output);
 
 	      if (verboseLevel > 3) {
 		G4cout << " After EquilibriumEvaporator " << G4endl;
-
 		output.printCollisionOutput();
 	      };
 
 	      TRFoutput.addOutgoingParticles(output.getOutgoingParticles());  
-	      TRFoutput.addTargetFragments(output.getNucleiFragments());         
+	      TRFoutput.addTargetFragments(output.getNucleiFragments());
 	    };
 	  };
 	 
-	  // convert to the LAB       
-	  G4bool withReflection = convertToTargetRestFrame.reflectionNeeded();       
-	  std::vector<G4InuclElementaryParticle> particles = 
-	    TRFoutput.getOutgoingParticles();
-
-	  if (!particles.empty()) { 
-	    particleIterator ipart;
-	    for(ipart = particles.begin(); ipart != particles.end(); ipart++) {
-	      G4CascadeMomentum mom = ipart->getMomentum();
-
-	      if (withReflection) mom[3] = -mom[3];
-	      mom = convertToTargetRestFrame.rotate(mom);
-	      ipart->setMomentum(mom); 
-	      mom = convertToTargetRestFrame.backToTheLab(ipart->getMomentum());
-	      ipart->setMomentum(mom); 
-	    };
-	    std::sort(particles.begin(), particles.end(), G4ParticleLargerEkin());
-	  };
-           
-	  std::vector<G4InuclNuclei> nucleus = TRFoutput.getNucleiFragments();
-
-	  if (!nucleus.empty()) { 
-	    nucleiIterator inuc;
-
-	    for (inuc = nucleus.begin(); inuc != nucleus.end(); inuc++) {
-	      G4CascadeMomentum mom = inuc->getMomentum(); 
-
-	      if (withReflection) mom[3] = -mom[3];
-	      mom = convertToTargetRestFrame.rotate(mom);
-	      inuc->setMomentum(mom);
-	      inuc->setEnergy(); 
-	      mom = convertToTargetRestFrame.backToTheLab(inuc->getMomentum());
-	      inuc->setMomentum(mom);
-	      inuc->setEnergy(); 
-	    };
-	  };
-	  globalOutput.addOutgoingParticles(particles);
-	  globalOutput.addTargetFragments(nucleus);
+	  // convert to the LAB
+	  TRFoutput.boostToLabFrame(convertToTargetRestFrame);
+
+	  globalOutput.addOutgoingParticles(TRFoutput.getOutgoingParticles());
+	  globalOutput.addTargetFragments(TRFoutput.getNucleiFragments());
 	  globalOutput.setOnShell(bullet, target);
-	  if(globalOutput.acceptable()) {
-
-	    return globalOutput;
-	  } else {
-	    globalOutput.reset();
-	  }; 
+	  if (globalOutput.acceptable()) return;
+
+	  globalOutput.reset();		// Clear and try again
 	};
 
@@ -243,149 +212,20 @@
 		 << itry_max << " attempts " << G4endl;
 	}
-
-	globalOutput.trivialise(bullet, target);
-
-	return globalOutput;        
       } else {
-
 	if (verboseLevel > 3) {
 	  G4cout << " InuclCollider -> inelastic interaction is impossible " << G4endl
 		 << " due to the coulomb barirer " << G4endl;
 	}
-
-	globalOutput.trivialise(bullet, target);
-
-	return globalOutput;
-      };
-	
+      }
+
+      globalOutput.trivialise(bullet, target);
+      return;
     } else {
-
       if (verboseLevel > 3) {
-	G4cout << " InuclCollider -> inter case " << intcase << G4endl;
+	G4cout << " InuclCollider -> inter case " << interCase.code() << G4endl;
       };
     };       
   };
 
-  return globalOutput;
+  return;
 }
-		     
-G4bool G4InuclCollider::inelasticInteractionPossible(G4InuclParticle* bullet,
-						     G4InuclParticle* target, 
-						     G4double ekin) const {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4InuclCollider::inelasticInteractionPossible" << G4endl;
-  }
-
-  const G4double coeff = 0.001 * 1.2;
-  const G4double one_third = 1.0 / 3.0;
-
-  G4bool possible = true;
-  G4double at;
-  G4double zt;
-  G4double ab;
-  G4double zb;
-
-  if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {
-    at = nuclei_target->getA();
-    zt = nuclei_target->getZ(); 
-    if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
-      ab = nuclei_bullet->getA();
-      zb = nuclei_bullet->getZ();     
-    } else {
-      G4InuclElementaryParticle* particle =
-	dynamic_cast<G4InuclElementaryParticle*>(bullet);
-
-      ab = 1;
-      zb = particle->getCharge();
-    }; 
-  } else {
-    if(G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) {
-      ab = nuclei_bullet->getA();
-      zb = nuclei_bullet->getZ();     
-
-      G4InuclElementaryParticle* particle =
-	dynamic_cast<G4InuclElementaryParticle*>(target);
-
-      at = 1;
-      zt = particle->getCharge();    
-    } else {
-
-      return possible;
-    };  
-  }; 
-
-  // VCOL used  for testing if elastic collision possible
-  G4double VCOL = coeff * zt * zb / (std::pow(at, one_third) + std::pow(ab, one_third)); 
-
-  // possible = VCOL < ekin; // NOTE: inelastic collision if not true
-  possible = true; // we force elastic
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4InuclCollider::inelasticInteractionPossible" << G4endl;
-    G4cout << " VCOL: " << VCOL << " ekin: " << ekin << " inelastic possible: " << possible << G4endl;
-  }
-
-  return possible;
-
-}
-	
-G4InteractionCase G4InuclCollider::bulletTargetSetter(G4InuclParticle* bullet,
-						      G4InuclParticle* target) const {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4InuclCollider::bulletTargetSetter" << G4endl;
-  }
-
-  G4InteractionCase interCase;
-
-  if (G4InuclNuclei* nuclei_target = dynamic_cast<G4InuclNuclei*>(target)) {     
-    if (G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet)) { // A + A         
-      interCase.setInterCase(2);
-      if (nuclei_target->getA() >= nuclei_bullet->getA()) {
-	interCase.setBulletTarget(bullet, target);
-      } else {
-	interCase.setBulletTarget(target, bullet);
-      }; 
-    } else {
-      interCase.setInterCase(1);
-      interCase.setBulletTarget(bullet, target);
-    }; 
-  } else {
-    G4InuclNuclei* nuclei_bullet = dynamic_cast<G4InuclNuclei*>(bullet);
-    if (nuclei_bullet) { 
-      G4InuclElementaryParticle* part =
-	dynamic_cast<G4InuclElementaryParticle*>(target);
-      if (part) {
-	interCase.setInterCase(1);
-	interCase.setBulletTarget(target, bullet);
-      };
-    }; 
-  };
-
-  return interCase;
-}       
-
-G4bool G4InuclCollider::explosion(G4InuclNuclei* target) const {
-
-  if (verboseLevel > 3) {
-    G4cout << " >>> G4InuclCollider::explosion" << G4endl;
-  }
-
-  const G4double a_cut = 20.0;
-  const G4double be_cut = 3.0;
-
-  G4double a = target->getA();
-  G4double z = target->getZ();
-  G4double eexs = target->getExitationEnergy();
-  G4bool explo = true;
-
-  if (a > a_cut) {
-    explo = false;
-  } else {
-    if (eexs < be_cut * bindingEnergy(a, z)) explo = false;
-  };   
-
-  return explo;
-}
- 
