Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc	(revision 1201)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc	(revision 1228)
@@ -787,11 +787,11 @@
 // get A and Z for the residual nucleus
   #if defined(debug_G4BinaryCascade) || defined(debug_BIC_GetExcitationEnergy)
-  G4int finalA = theTargetList.size()+theCapturedList.size();
-  G4int finalZ = GetTotalCharge(theTargetList)+GetTotalCharge(theCapturedList);
-  if ( (currentA - finalA) != 0 || (currentZ - finalZ) != 0 )
-  {
-     G4cerr << "G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} " 
-            << currentA << " " << finalA << " "<< currentZ << " " << finalZ << G4endl;
-  }
+     G4int finalA = theTargetList.size()+theCapturedList.size();
+     G4int finalZ = GetTotalCharge(theTargetList)+GetTotalCharge(theCapturedList);
+     if ( (currentA - finalA) != 0 || (currentZ - finalZ) != 0 )
+     {
+	G4cerr << "G4BIC:GetExcitationEnergy(): Nucleon counting error current/final{A,Z} " 
+               << currentA << " " << finalA << " "<< currentZ << " " << finalZ << G4endl;
+     }
   
   #endif
@@ -852,5 +852,5 @@
       initialExc = theInitial4Mom.mag()-
            G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(Z, A);
-	   G4cout << " Initial nucleus A Z" << A << " " << Z << initialExc << G4endl; 
+	   G4cout << "GetExcitationEnergy: Initial nucleus A Z " << A << " " << Z << " " << initialExc << G4endl; 
     }
   }
@@ -1081,5 +1081,4 @@
   G4int initialBaryon(0);
   G4int initialCharge(0);
-  G4double initial_Efermi(0);
   G4LorentzVector mom4Primary(0);
   
@@ -1090,24 +1089,27 @@
   }
 
-// for primary resonances, subtract neutron ( = proton) field ( ie. add std::abs(field))
   G4int PDGcode=std::abs(primary->GetDefinition()->GetPDGEncoding());
   mom4Primary=primary->Get4Momentum();
-  initial_Efermi=RKprop->GetField(primary->GetDefinition()->GetPDGEncoding(),primary->GetPosition());
-
-  if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
-  {
-     initial_Efermi = RKprop->GetField(G4Neutron::Neutron()->GetPDGEncoding(),
-                                                primary->GetPosition());
-     primary->Update4Momentum(mom4Primary.e() - initial_Efermi);
-  }
-
-  std::vector<G4KineticTrack *>::iterator titer;
-  for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
-  {
-     G4ParticleDefinition * aDef=(*titer)->GetDefinition();
-     G4int aCode=aDef->GetPDGEncoding();
-     G4ThreeVector aPos=(*titer)->GetPosition();
-     initial_Efermi+= RKprop->GetField(aCode, aPos);
-//     initial_Efermi+= RKprop->GetField((*titer)->GetDefinition()->GetPDGEncoding(),(*titer)->GetPosition());
+
+// for primary resonances, subtract neutron ( = proton) field ( ie. add std::abs(field))
+  G4double initial_Efermi(0);
+  if (primary->GetState() == G4KineticTrack::inside ) {
+     initial_Efermi=RKprop->GetField(primary->GetDefinition()->GetPDGEncoding(),primary->GetPosition());
+
+     if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
+     {
+	initial_Efermi = RKprop->GetField(G4Neutron::Neutron()->GetPDGEncoding(),
+                                                   primary->GetPosition());
+	primary->Update4Momentum(mom4Primary.e() - initial_Efermi);
+     }
+
+     std::vector<G4KineticTrack *>::iterator titer;
+     for ( titer=target_collection.begin() ; titer!=target_collection.end(); ++titer)
+     {
+	G4ParticleDefinition * aDef=(*titer)->GetDefinition();
+	G4int aCode=aDef->GetPDGEncoding();
+	G4ThreeVector aPos=(*titer)->GetPosition();
+	initial_Efermi+= RKprop->GetField(aCode, aPos);
+     }
   }
 //****************************************
@@ -1140,5 +1142,5 @@
      G4bool lateParticleCollision= (!haveTarget) && products && products->size() == 1;
 	//  if ( lateParticleCollision ) G4cout << " Added late particle--------------------------"<<G4endl;
-	//  if ( products ) PrintKTVector(products, " reaction products");
+	//  if ( lateParticleCollision && products ) PrintKTVector(products, " reaction products");
 //****************************************  
 
@@ -1172,38 +1174,40 @@
   }
 
-   G4double final_Efermi(0);
-   G4KineticTrackVector resonances;
-   for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
-   {
-       G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
-//       G4cout << " PDGcode, state " << PDGcode << " " << (*i)->GetState()<<G4endl;
-       final_Efermi+=RKprop->GetField(PDGcode,(*i)->GetPosition());
-       if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
-       {  
-	  resonances.push_back(*i);
-       }
-   }	
-   if ( resonances.size() > 0 ) 
-   {  
-      G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
-      for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
-      {
-	  G4LorentzVector mom=(*res)->Get4Momentum();
-	  G4double mass2=mom.mag2();
-	  G4double newEnergy=mom.e() + delta_Fermi;
-	  G4double newEnergy2= newEnergy*newEnergy;
-	  	//G4cout << "mom = " << mom <<" newE " << newEnergy<< G4endl;
-	  if ( newEnergy2 < mass2 )
-	  {
-             ClearAndDestroy(products);
-             if (target_collection.size() == 0 ) FindDecayCollision(primary);  // for decay, sample new decay
-	     delete products;
-	     return false;
-	  }
-//	  G4cout << " correct resonance from /to " << mom.e() << " / " << newEnergy<< G4endl;
-	  G4ThreeVector mom3=std::sqrt(newEnergy2 - mass2) * mom.vect().unit();
-	  (*res)->Set4Momentum(G4LorentzVector(mom3,newEnergy));
-      }
-   }
+  if (primary->GetState() == G4KineticTrack::inside ) {   // if the primary was outside, nothing to correct
+     G4double final_Efermi(0);
+     G4KineticTrackVector resonances;
+     for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
+     {
+	 G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
+         //       G4cout << " PDGcode, state " << PDGcode << " " << (*i)->GetState()<<G4endl;
+	 final_Efermi+=RKprop->GetField(PDGcode,(*i)->GetPosition());
+	 if ( PDGcode > 1000 && PDGcode != 2112 && PDGcode != 2212 )
+	 {  
+	    resonances.push_back(*i);
+	 }
+     }	
+     if ( resonances.size() > 0 ) 
+     {  
+	G4double delta_Fermi= (initial_Efermi-final_Efermi)/resonances.size();
+	for (std::vector<G4KineticTrack *>::iterator res=resonances.begin(); res != resonances.end(); res++)
+	{
+	    G4LorentzVector mom=(*res)->Get4Momentum();
+	    G4double mass2=mom.mag2();
+	    G4double newEnergy=mom.e() + delta_Fermi;
+	    G4double newEnergy2= newEnergy*newEnergy;
+	  	  //G4cout << "mom = " << mom <<" newE " << newEnergy<< G4endl;
+	    if ( newEnergy2 < mass2 )
+	    {
+               ClearAndDestroy(products);
+               if (target_collection.size() == 0 ) FindDecayCollision(primary);  // for decay, sample new decay
+	       delete products;
+	       return false;
+	    }
+            //	  G4cout << " correct resonance from /to " << mom.e() << " / " << newEnergy<< G4endl;
+	    G4ThreeVector mom3=std::sqrt(newEnergy2 - mass2) * mom.vect().unit();
+	    (*res)->Set4Momentum(G4LorentzVector(mom3,newEnergy));
+	}
+     }
+  }
 
 #ifdef debug_BIC_ApplyCollision
@@ -1218,7 +1222,9 @@
     if ( ! lateParticleCollision ) 
     {
-       (*i)->SetState(primary->GetState());  // secondaries are created inside nucleus, except for decay in propagate
-       finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
-       finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge());
+       (*i)->SetState(primary->GetState());  // decay may be anywhere!
+       if ( (*i)->GetState() == G4KineticTrack::inside ){
+          finalBaryon+=(*i)->GetDefinition()->GetBaryonNumber();
+          finalCharge+=G4lrint((*i)->GetDefinition()->GetPDGCharge());
+       }
     } else {
        G4double tin=0., tout=0.; 
@@ -1247,5 +1253,5 @@
        }
        
-//       G4cout << " PDGcode, state " << (*i)->GetDefinition()->GetPDGEncoding() << " " << (*i)->GetState()<<G4endl;
+       //G4cout << " PDGcode, state " << (*i)->GetDefinition()->GetPDGEncoding() << " " << (*i)->GetState()<<G4endl;
 
     }   
@@ -1269,5 +1275,5 @@
   currentA += finalBaryon-initialBaryon;
   currentZ += finalCharge-initialCharge;
-	//G4cout << " currentA, Z aft: " << currentA << " " << currentZ << G4endl;
+	//G4cout << " ApplyCollision currentA, Z aft: " << currentA << " " << currentZ << G4endl;
   
   G4KineticTrackVector oldSecondaries;
@@ -2351,5 +2357,5 @@
      G4cerr << "G4BinaryCascade::GetFinal4Momentum - Fatal"<<G4endl;
      G4KineticTrackVector::iterator i;
-     G4cerr <<" Initial nucleus "<<theInitial4Mom<<G4endl;
+     G4cerr <<" GetFinal4Momentum: Initial nucleus "<<theInitial4Mom<<G4endl;
      for(i = theProjectileList.begin() ; i != theProjectileList.end(); ++i)
      {
