Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4Absorber.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4Absorber.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4Absorber.cc	(revision 1196)
@@ -56,6 +56,6 @@
   if (kt.Get4Momentum().e() - kt.GetActualMass() < theCutOnP) 
   {
-      if(kt.GetDefinition() == G4PionPlus::PionPlus() ||
-	 kt.GetDefinition() == G4PionZero::PionZero() ||
+      if(kt.GetDefinition() == G4PionPlus::PionPlus()  ||
+	 kt.GetDefinition() == G4PionZero::PionZero()  ||
 	 kt.GetDefinition() == G4PionMinus::PionMinus())
       {
Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc	(revision 1196)
@@ -97,9 +97,9 @@
   theLateParticle= new G4BCLateParticle;
   theImR.push_back(theDecay);
+  G4MesonAbsorption * aAb=new G4MesonAbsorption;
+  theImR.push_back(aAb);
   G4Scatterer * aSc=new G4Scatterer;
   theH1Scatterer = new G4Scatterer;
   theImR.push_back(aSc);
-  G4MesonAbsorption * aAb=new G4MesonAbsorption;
-  theImR.push_back(aAb);
 
   thePropagator = new G4RKPropagation;
@@ -107,5 +107,5 @@
   theBCminP = 45*MeV;
   theCutOnP = 90*MeV;
-  theCutOnPAbsorb= 0*MeV;
+  theCutOnPAbsorb= 0*MeV;   // No Absorption of slow Mesons, other than above G4MesonAbsorption
 
   theExcitationHandler = new G4ExcitationHandler;
@@ -155,4 +155,5 @@
   eventcounter++;
   if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number starts ######### "<<eventcounter<<G4endl;
+
   G4LorentzVector initial4Momentum = aTrack.Get4Momentum();
   G4ParticleDefinition * definition = const_cast<G4ParticleDefinition *>(aTrack.GetDefinition());
@@ -238,9 +239,4 @@
   theParticleChange.SetStatusChange(stopAndKill);
   G4ReactionProductVector::iterator iter;
-  G4double Efinal=0;
-  	if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
-	{
-	   G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl;
-	}
 
   for(iter = products->begin(); iter != products->end(); ++iter)
@@ -253,16 +249,8 @@
     theParticleChange.AddSecondary(aNew);
 
-//   G4cout << " Secondary E - Ekin / p " <<
-//      (*iter)->GetDefinition()->GetParticleName() << " " <<
-//      (*iter)->GetTotalEnergy() << " - " <<
-//      (*iter)->GetKineticEnergy()<< " / " <<
-//      (*iter)->GetMomentum().x() << " " <<
-//      (*iter)->GetMomentum().y() << " " <<
-//      (*iter)->GetMomentum().z() << G4endl;
-      Efinal += (*iter)->GetTotalEnergy();
-  }
-
-//  G4cout << "e outgoing/ total : " << Efinal << " " << Efinal+GetFinal4Momentum().e()<< G4endl;
-
+  }
+
+  // DebugEpConservation(track, products);
+              
   ClearAndDestroy(products);
   delete products;
@@ -325,8 +313,17 @@
   if(nucleus->GetMass()>120) theCutOnP = 45*MeV;
 
+  G4double StartingTime=DBL_MAX;        // Search for minimal formation time Uzhi
+  for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)    // Uzhi
+  {                                                                       // Uzhi
+    if((*iter)->GetFormationTime() < StartingTime)                        // Uzhi
+        StartingTime = (*iter)->GetFormationTime();                       // Uzhi
+  }                                                                       // Uzhi
+
   for(iter = secondaries->begin(); iter != secondaries->end(); ++iter)
   {
 //  G4cout << " Formation time : " << (*iter)->GetDefinition()->GetParticleName() << " " 
 // 	 << (*iter)->GetFormationTime() << G4endl;
+    G4double FormTime = (*iter)->GetFormationTime() - StartingTime;       // Uzhi
+    (*iter)->SetFormationTime(FormTime);                                  // Uzhi
     if( (*iter)->GetState() == G4KineticTrack::undefined ) 
     {
@@ -363,5 +360,5 @@
   G4bool haveProducts = false;
   G4int collisionCount=0;
-  theMomentumTransfer=0;
+  theMomentumTransfer=G4ThreeVector(0,0,0);
   while(theCollisionMgr->Entries() > 0)
   {
@@ -394,5 +391,5 @@
        if (!DoTimeStep(nextCollision->GetCollisionTime()-theCurrentTime) )
        {
-	   // Check if nextCollision is still valid, ie. partcile did not leave nucleus
+	   // Check if nextCollision is still valid, ie. particle did not leave nucleus
 	   if (theCollisionMgr->GetNextCollision() != nextCollision )
 	   {
@@ -539,81 +536,142 @@
   }
 
-
 // find a fragment and call the precompound model.
   G4Fragment * fragment = 0;
   G4ReactionProductVector * precompoundProducts = 0;
- 
-   if ( ExcitationEnergy >= 0 ) 
-   {
-       fragment = FindFragments();
-
-    //  theDeExcitation =0;
-       if(fragment && fragment->GetA() >1.5) // hpw @@@@ still to add the nucleon, if one exists.
+
+  G4LorentzVector pFragment(0);
+        // G4cout << " final4mon " << GetFinal4Momentum() /MeV << G4endl; 
+
+//   if ( ExcitationEnergy >= 0 )                                         // closed by Uzhi
+//   {                                                                    // closed by Uzhi
+  fragment = FindFragments();
+  if(fragment)                                                            // Uzhi
+  {                                                                       // Uzhi
+       if(fragment->GetA() >1.5)                                          // Uzhi
        {
-	 if (theDeExcitation)
+	 if (theDeExcitation)                // pre-compound
 	 {
+	      // G4cout << " going to preco with fragment 4 mom " << fragment->GetMomentum() << G4endl;
+             pFragment=fragment->GetMomentum();
              precompoundProducts= theDeExcitation->DeExcite(*fragment);
              delete fragment;
              fragment=0;
-	 } else if (theExcitationHandler)
+	 } else if (theExcitationHandler)    // de-excitation
 	 {
+	      // G4cout << " going to de-excit with fragment 4 mom " << fragment->GetMomentum() << G4endl;
+             pFragment=fragment->GetMomentum();
 	     precompoundProducts=theExcitationHandler->BreakItUp(*fragment);
              delete fragment;
              fragment=0;
 	 }
-       } else {
+       } else 
+       {                                   // fragment->GetA() < 1.5
+	 precompoundProducts = new G4ReactionProductVector();
+	 std::vector<G4KineticTrack *>::iterator i;
          if ( theTargetList.size() == 1 )
 	 {
-	     precompoundProducts = new G4ReactionProductVector();
-	     std::vector<G4KineticTrack *>::iterator i=theTargetList.begin();
+             i=theTargetList.begin();
 	     G4ReactionProduct * aNew = new G4ReactionProduct((*i)->GetDefinition());
 	     aNew->SetTotalEnergy((*i)->GetDefinition()->GetPDGMass());       
-	     aNew->SetMomentum(G4ThreeVector(0));		// see boost fro preCompoundProducts below..
+	     aNew->SetMomentum(G4ThreeVector(0));// see boost for preCompoundProducts below..
 	     precompoundProducts->push_back(aNew);
-	  } else if ( theTargetList.size() > 1)
+	 } 
+
+         if ( theCapturedList.size() == 1 )                               // Uzhi
+	 {                                                                // Uzhi
+             i=theCapturedList.begin();                                   // Uzhi
+	     G4ReactionProduct * aNew = new G4ReactionProduct((*i)->GetDefinition()); // Uzhi
+	     aNew->SetTotalEnergy((*i)->GetDefinition()->GetPDGMass());   // Uzhi
+	     aNew->SetMomentum(G4ThreeVector(0));// see boost below..     // Uzhi
+	     precompoundProducts->push_back(aNew);                        // Uzhi
+	 }                                                                // Uzhi
+       }                            // End of fragment->GetA() < 1.5
+  } else                            // End of if(fragment)
+  {                                 // No fragment, can be neutrons only  // Uzhi
+     precompoundProducts = new G4ReactionProductVector();                      
+     
+     if ( (theTargetList.size()+theCapturedList.size()) > 0 )
+     {                                                                      
+	std::vector<G4KineticTrack *>::iterator aNuc;                             
+	G4LorentzVector aVec;                                                     
+	std::vector<G4double> masses;                                             
+	G4double sumMass(0);
+
+	if ( theTargetList.size() != 0)                                      // Uzhi
+	{
+           for ( aNuc=theTargetList.begin(); aNuc != theTargetList.end(); aNuc++)
+           {
+              G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();
+              masses.push_back(mass);
+              sumMass += mass;
+           }
+	}                                                                    // Uzhi
+
+	if ( theCapturedList.size() != 0)                                    // Uzhi
+	{                                                                    // Uzhi
+           for(aNuc = theCapturedList.begin();                               // Uzhi
+               aNuc != theCapturedList.end(); aNuc++)                        // Uzhi
+           {                                                                 // Uzhi
+              G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();          // Uzhi
+              masses.push_back(mass);                                        // Uzhi
+              sumMass += mass;                                               // Uzhi
+           }                 
+	}
+
+	G4LorentzVector finalP=GetFinal4Momentum();
+	G4FermiPhaseSpaceDecay decay;
+	// G4cout << " some neutrons? " << masses.size() <<" " ;
+	// G4cout<< theTargetList.size()<<" "<<finalP <<" " << finalP.mag()<<G4endl;
+
+	G4double eCMS=finalP.mag();
+	if ( eCMS < sumMass )                    // @@GF --- Cheat!!
+	{
+           eCMS=sumMass + (2*MeV*masses.size());     
+	   finalP.setE(std::sqrt(finalP.vect().mag2() + sqr(eCMS)));
+	}
+
+	precompoundLorentzboost.set(finalP.boostVector());
+	std::vector<G4LorentzVector*> * momenta=decay.Decay(eCMS,masses);
+	std::vector<G4LorentzVector*>::iterator aMom=momenta->begin();
+
+	if ( theTargetList.size() != 0)
+	{
+	  for ( aNuc=theTargetList.begin(); 
+               (aNuc != theTargetList.end()) && (aMom!=momenta->end()); 
+        	aNuc++, aMom++ )
 	  {
-	     precompoundProducts = new G4ReactionProductVector();
-	     std::vector<G4KineticTrack *>::iterator aNuc;
-	     G4LorentzVector aVec;
-	     std::vector<G4double> masses;
-	     G4double sumMass(0);
-	     for ( aNuc=theTargetList.begin(); aNuc != theTargetList.end(); aNuc++)
-	     {
-	        G4double mass=(*aNuc)->GetDefinition()->GetPDGMass();
-		masses.push_back(mass);
-		sumMass += mass;
-	     }
-	     G4LorentzVector finalP=GetFinal4Momentum();
-	     G4FermiPhaseSpaceDecay decay;
-//	     G4cout << " some neutrons? " << masses.size() <<" " << theTargetList.size()<<" "<<finalP <<" " << finalP.mag()<<G4endl;
-             G4double eCMS=finalP.mag();
-	     if ( eCMS < sumMass )                    // @@GF --- Cheat!!
-	     {
-	        eCMS=sumMass + (2*MeV*masses.size());     
-		finalP.setE(std::sqrt(finalP.vect().mag2() + sqr(eCMS)));
-	     }
-	     precompoundLorentzboost.set(finalP.boostVector());
-	     std::vector<G4LorentzVector*> * momenta=decay.Decay(eCMS,masses);
-	     std::vector<G4LorentzVector*>::iterator aMom=momenta->begin();
-	     for ( aNuc=theTargetList.begin(); 
-	          (aNuc != theTargetList.end()) && (aMom!=momenta->end()); 
-		  aNuc++, aMom++ )
-	     {
-	        G4ReactionProduct * aNew = new G4ReactionProduct((*aNuc)->GetDefinition());
-		aNew->SetTotalEnergy((*aMom)->e());
-		aNew->SetMomentum((*aMom)->vect());
-		precompoundProducts->push_back(aNew);
-//		G4cout << " only neutrons? " <<  (*aNuc)->GetDefinition()->GetParticleName() << G4endl;
-                delete *aMom;
-	     }
-	     delete momenta;
+             G4ReactionProduct * aNew = new G4ReactionProduct((*aNuc)->GetDefinition());
+             aNew->SetTotalEnergy((*aMom)->e());
+             aNew->SetMomentum((*aMom)->vect());
+             precompoundProducts->push_back(aNew);
+
+             delete *aMom;
 	  }
-       }
-
-   } 
+	}
+
+	if ( theCapturedList.size() != 0)                                    // Uzhi
+	{                                                                    // Uzhi
+	  for ( aNuc=theCapturedList.begin();                                // Uzhi
+               (aNuc != theCapturedList.end()) && (aMom!=momenta->end());    // Uzhi
+		aNuc++, aMom++ )                                             // Uzhi
+	  {                                                                  // Uzhi
+             G4ReactionProduct * aNew = new G4ReactionProduct(               // Uzhi
+                                       (*aNuc)->GetDefinition());            // Uzhi
+             aNew->SetTotalEnergy((*aMom)->e());                             // Uzhi
+             aNew->SetMomentum((*aMom)->vect());                             // Uzhi
+             precompoundProducts->push_back(aNew);                           // Uzhi
+             delete *aMom;                                                   // Uzhi
+	  }                                                                  // Uzhi
+	}                                                                    // Uzhi
+
+	if(momenta) delete momenta;
+     }  
+  }                   // End if(!fragment)
+
 
   {
 // fill in products the outgoing particles
      G4double Ekinout=0;
+     G4LorentzVector pSumBic(0);
      size_t i(0);
      for(i = 0; i< theFinalState.size(); i++)
@@ -636,4 +694,5 @@
 	 aNew->SetTotalEnergy(kt->Get4Momentum().e());
 	 Ekinout += aNew->GetKineticEnergy();
+	 pSumBic += kt->Get4Momentum();
 	 if(kt->IsParticipant()) 
 	 {
@@ -655,5 +714,5 @@
 	     } else
 	     {
-		G4cout << "final un stable : ";
+		G4cout << "final non stable : ";
 	     }
 	     G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl;
@@ -666,4 +725,5 @@
   }
 // add precompound products to products
+  G4LorentzVector pSumPreco(0), pPreco(0);
   if ( precompoundProducts )
   {
@@ -673,4 +733,5 @@
 // boost back to system of moving nucleus
          G4LorentzVector pProduct((*j)->GetMomentum(),(*j)->GetTotalEnergy());
+	 pPreco+= pProduct;
 #ifdef debug_BIC_Propagate_finals
 	 G4cout << " pProduct be4 boost " <<pProduct << G4endl;
@@ -680,4 +741,5 @@
 	 G4cout << " pProduct aft boost " <<pProduct << G4endl;
 #endif
+         pSumPreco += pProduct;
          (*j)->SetTotalEnergy(pProduct.e());
          (*j)->SetMomentum(pProduct.vect());
@@ -685,5 +747,6 @@
 	 products->push_back(*j);
        }
-
+	 // G4cout << " unboosted preco result mom " << pPreco / MeV << "  ..- fragmentMom " << (pPreco - pFragment)/MeV<< G4endl;
+	 // G4cout << " preco result mom " << pSumPreco / MeV << "  ..-file4Mom " << (pSumPreco - GetFinal4Momentum())/MeV<< G4endl;
        precompoundProducts->clear();
        delete precompoundProducts;
@@ -740,8 +803,10 @@
      nucleusMass = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA);
   } 
-  else if (currentZ==0 && currentA==1 )
-  {
-     nucleusMass = G4Neutron::Neutron()->GetPDGMass();
-  } 
+  else if (currentZ==0 )     // Uzhi && currentA==1 )                     // Uzhi
+  {                                                                       // Uzhi
+     if(currentA == 1) {nucleusMass = G4Neutron::Neutron()->GetPDGMass();}// Uzhi
+     else              {nucleusMass = GetFinalNucleusMomentum().mag()     // Uzhi
+                                      - 3.*MeV*currentA;}                 // Uzhi
+  }                                                                       // Uzhi
   else
   {
@@ -929,4 +994,5 @@
 //----------------------------------------------------------------------------
 {
+    G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
     if ( secondary->GetTrackingMomentum().mag2() < -1.*eV )
     {
@@ -952,4 +1018,13 @@
     }
 
+    //  for barions, correct for fermi energy
+    G4int PDGcode=std::abs(secondary->GetDefinition()->GetPDGEncoding());
+    if ( PDGcode > 1000 ) 
+    {
+       G4double initial_Efermi = RKprop->GetField(G4Neutron::Neutron()->GetPDGEncoding(),
+                                          secondary->GetPosition()); 
+       secondary->Update4Momentum(secondary->Get4Momentum().e() - initial_Efermi);
+    }
+
 #ifdef debug_BIC_FindCollision
     G4cout << "FindLateP Particle, 4-mom, times newState " 
@@ -1019,9 +1094,12 @@
   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());
+     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)
@@ -1039,7 +1117,13 @@
 
   #ifdef debug_BIC_ApplyCollision
-      if ( !products )
+        G4bool havePion=false;
+        for ( std::vector<G4KineticTrack *>::iterator i =products->begin(); i != products->end(); i++)
+   	{
+       		G4int PDGcode=std::abs((*i)->GetDefinition()->GetPDGEncoding());
+		if (std::abs(PDGcode)==211 || PDGcode==111 ) havePion=true;
+	}	
+     if ( !products  || havePion)
       {
-            G4cout << " Collision type: "<< typeid(*collision->GetGenerator()).name()
+            G4cout << " Collision " << collision << ", type: "<< typeid(*collision->GetGenerator()).name()
 		<< ", with NO products! " <<G4endl;
 	   G4cout << G4endl<<"Initial condition are these:"<<G4endl;
@@ -1123,6 +1207,6 @@
    }
 
-#ifdef debug_BIC_ApplyCollisionXX
-  DebugApplyCollsion(collision, products);
+#ifdef debug_BIC_ApplyCollision
+  DebugApplyCollision(collision, products);
 #endif
 
@@ -1231,96 +1315,4 @@
 }
 
-//----------------------------------------------------------------------------
-void G4BinaryCascade::DebugApplyCollision(G4CollisionInitialState * collision, 
-                                          G4KineticTrackVector * products)
-{
-//**-------------------------begin debug Block---------------------------
-
-  G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
-
-  G4KineticTrackVector debug1;
-  debug1.push_back(collision->GetPrimary());
-  PrintKTVector(&debug1,std::string(" Primary particle"));
-  PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
-  PrintKTVector(products,std::string(" Scatterer products"));
-  
-  G4double thisExcitation(0);
-//  excitation energy from this collision
-//  initial state:
-  G4double initial(0);
-  G4KineticTrack * kt=collision->GetPrimary();
-  initial +=  kt->Get4Momentum().e();
-  
-  initial +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
-  initial -=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
-  G4cout << "prim. E/field/Barr/Sum " << kt->Get4Momentum().e()
-          << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
-          << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 
-	  << " " << initial << G4endl;;
-  
-  G4KineticTrackVector ktv=collision->GetTargetCollection();
-  for ( unsigned int it=0; it < ktv.size(); it++)
-  {
-     kt=ktv[it];
-     initial +=  kt->Get4Momentum().e();
-     thisExcitation += kt->GetDefinition()->GetPDGMass() 
-     		     - kt->Get4Momentum().e() 
-		     - RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
-//     initial +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
-//     initial -=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
-  G4cout << "Targ. def/E/field/Barr/Sum " <<  kt->GetDefinition()->GetPDGEncoding()
-  	  << " " << kt->Get4Momentum().e()
-          << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
-          << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 
-	  << " " << initial <<" Excit " << thisExcitation << G4endl;;
-  }
-  
-  G4double final(0);
-  G4double mass_out(0);
-  G4int product_barions(0);
-  if ( products ) 
-  {
-     for ( unsigned int it=0; it < products->size(); it++)
-     {
-	kt=(*products)[it];
-	final +=  kt->Get4Momentum().e();
-	final +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
-	final +=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
-	if ( kt->GetDefinition()->GetBaryonNumber()==1 ) product_barions++;
-	mass_out += kt->GetDefinition()->GetPDGMass();
-     G4cout << "sec. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
-  	     << " " << kt->Get4Momentum().e()
-             << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
-             << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 
-	     << " " << final << G4endl;;
-     }
-  }
-
-
-  G4int finalA = currentA;
-  G4int finalZ = currentZ;
-  if ( products )
-  {
-     finalA -= product_barions;
-     finalZ -= GetTotalCharge(*products);
-  }
-  G4double delta = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA) 
-                   - (G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA) 
-		       + mass_out); 
-  G4cout << " current/final a,z " << currentA << " " << currentZ << " "<< finalA<< " "<< finalZ 
-        <<  " delta-mass " << delta<<G4endl;
-  final+=delta;
-    mass_out  = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA);
-  G4cout << " initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy
-         << " " <<   final << " "
-	 <<  mass_out<<" " 
-	 <<  currentInitialEnergy - final - mass_out
-	 << G4endl;
-   currentInitialEnergy-=final;	 
-
-//**-------------------------end debug Block---------------------------
-
-}
-//----------------------------------------------------------------------------
 
 //----------------------------------------------------------------------------
@@ -2706,9 +2698,124 @@
 }
 
-
-
-
-
-
-
-
+//----------------------------------------------------------------------------
+void G4BinaryCascade::DebugApplyCollision(G4CollisionInitialState * collision, 
+                                          G4KineticTrackVector * products)
+{
+  G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
+
+  G4KineticTrackVector debug1;
+  debug1.push_back(collision->GetPrimary());
+  PrintKTVector(&debug1,std::string(" Primary particle"));
+  PrintKTVector(&collision->GetTargetCollection(),std::string(" Target particles"));
+  PrintKTVector(products,std::string(" Scatterer products"));
+  
+  G4double thisExcitation(0);
+//  excitation energy from this collision
+//  initial state:
+  G4double initial(0);
+  G4KineticTrack * kt=collision->GetPrimary();
+  initial +=  kt->Get4Momentum().e();
+  
+  initial +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
+  initial -=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
+  G4cout << "prim. E/field/Barr/Sum " << kt->Get4Momentum().e()
+          << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
+          << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 
+	  << " " << initial << G4endl;;
+  
+  G4KineticTrackVector ktv=collision->GetTargetCollection();
+  for ( unsigned int it=0; it < ktv.size(); it++)
+  {
+     kt=ktv[it];
+     initial +=  kt->Get4Momentum().e();
+     thisExcitation += kt->GetDefinition()->GetPDGMass() 
+     		     - kt->Get4Momentum().e() 
+		     - RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
+//     initial +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
+//     initial -=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
+  G4cout << "Targ. def/E/field/Barr/Sum " <<  kt->GetDefinition()->GetPDGEncoding()
+  	  << " " << kt->Get4Momentum().e()
+          << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
+          << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 
+	  << " " << initial <<" Excit " << thisExcitation << G4endl;;
+  }
+  
+  G4double final(0);
+  G4double mass_out(0);
+  G4int product_barions(0);
+  if ( products ) 
+  {
+     for ( unsigned int it=0; it < products->size(); it++)
+     {
+	kt=(*products)[it];
+	final +=  kt->Get4Momentum().e();
+	final +=  RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition());
+	final +=  RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding());
+	if ( kt->GetDefinition()->GetBaryonNumber()==1 ) product_barions++;
+	mass_out += kt->GetDefinition()->GetPDGMass();
+     G4cout << "sec. def/E/field/Barr/Sum " << kt->GetDefinition()->GetPDGEncoding()
+  	     << " " << kt->Get4Momentum().e()
+             << " " << RKprop->GetField(kt->GetDefinition()->GetPDGEncoding(),kt->GetPosition())
+             << " " << RKprop->GetBarrier(kt->GetDefinition()->GetPDGEncoding()) 
+	     << " " << final << G4endl;;
+     }
+  }
+
+
+  G4int finalA = currentA;
+  G4int finalZ = currentZ;
+  if ( products )
+  {
+     finalA -= product_barions;
+     finalZ -= GetTotalCharge(*products);
+  }
+  G4double delta = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(currentZ,currentA) 
+                   - (G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA) 
+		       + mass_out); 
+  G4cout << " current/final a,z " << currentA << " " << currentZ << " "<< finalA<< " "<< finalZ 
+        <<  " delta-mass " << delta<<G4endl;
+  final+=delta;
+    mass_out  = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(finalZ,finalA);
+  G4cout << " initE/ E_out/ Mfinal/ Excit " << currentInitialEnergy
+         << " " <<   final << " "
+	 <<  mass_out<<" " 
+	 <<  currentInitialEnergy - final - mass_out
+	 << G4endl;
+   currentInitialEnergy-=final;	 
+
+}
+
+//----------------------------------------------------------------------------
+void G4BinaryCascade::DebugEpConservation(const G4HadProjectile & aTrack,
+					  G4ReactionProductVector* products)			   
+{  
+  G4ReactionProductVector::iterator iter;
+  G4double Efinal(0);
+  G4ThreeVector pFinal(0);
+  	if (std::abs(theParticleChange.GetWeightChange() -1 ) > 1e-5 )
+	{
+	   G4cout <<" BIC-weight change " << theParticleChange.GetWeightChange()<< G4endl;
+	}
+
+  for(iter = products->begin(); iter != products->end(); ++iter)
+  {
+
+//   G4cout << " Secondary E - Ekin / p " <<
+//      (*iter)->GetDefinition()->GetParticleName() << " " <<
+//      (*iter)->GetTotalEnergy() << " - " <<
+//      (*iter)->GetKineticEnergy()<< " / " <<
+//      (*iter)->GetMomentum().x() << " " <<
+//      (*iter)->GetMomentum().y() << " " <<
+//      (*iter)->GetMomentum().z() << G4endl;
+      Efinal += (*iter)->GetTotalEnergy();
+      pFinal += (*iter)->GetMomentum();
+  }
+
+//  G4cout << "e outgoing/ total : " << Efinal << " " << Efinal+GetFinal4Momentum().e()<< G4endl;
+    G4cout << "BIC E/p delta " << 
+    (aTrack.Get4Momentum().e()+ the3DNucleus->GetMass() - Efinal)/MeV <<
+    " MeV / mom " << (aTrack.Get4Momentum()  - pFinal ) /MeV << G4endl;
+
+}
+
+//----------------------------------------------------------------------------
Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryLightIonReaction.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryLightIonReaction.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryLightIonReaction.cc	(revision 1196)
@@ -38,11 +38,15 @@
 #include <cmath>
  
-  G4BinaryLightIonReaction::G4BinaryLightIonReaction()
-    : G4HadronicInteraction("Binary Cascade"), theModel(), theHandler(), 
-      theProjectileFragmentation(&theHandler) {}
+G4BinaryLightIonReaction::G4BinaryLightIonReaction()
+    : G4HadronicInteraction("Binary Cascade"), theModel() , 
+      theHandler(0) , theProjectileFragmentation(0)
+{
+    theHandler= new G4ExcitationHandler; 
+    SetPrecompound(new G4PreCompoundModel(theHandler));
+}
   
-  G4HadFinalState *G4BinaryLightIonReaction::
+G4HadFinalState *G4BinaryLightIonReaction::
   ApplyYourself(const G4HadProjectile &aTrack, G4Nucleus & targetNucleus )
-  { 
+{ 
   static G4int eventcounter=0;
   eventcounter++;
@@ -118,5 +122,5 @@
 //      G4cout << "Fragment INFO "<< a1+a2 <<" "<<z1+z2<<" "
 //             << aL <<" "<<preFragDef->GetParticleName()<<G4endl;
-      cascaders = theProjectileFragmentation.DeExcite(aPreFrag);
+      cascaders = theProjectileFragmentation->DeExcite(aPreFrag);
       G4double tSum = 0;
       for(size_t count = 0; count<cascaders->size(); count++)
@@ -402,5 +406,5 @@
 	resDef = G4ParticleTable::GetParticleTable()->FindIon(resZ,resA,0,resZ);  
 	aProRes.SetParticleDefinition(resDef);
-	proFrag = theHandler.BreakItUp(aProRes);
+	proFrag = theHandler->BreakItUp(aProRes);
       if ( momentum.vect().mag() > momentum.e() )
       {
Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4GeneratorPrecompoundInterface.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4GeneratorPrecompoundInterface.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4GeneratorPrecompoundInterface.cc	(revision 1196)
@@ -31,4 +31,8 @@
 // HPW, 10DEC 98, the decay part originally written by Gunter Folger in his FTF-test-program.
 //
+   
+    G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface()
+    : CaptureThreshold(80*MeV)
+    {}
       
    G4HadFinalState* G4GeneratorPrecompoundInterface::
@@ -51,5 +55,5 @@
      result1=theSecondaries;
      result=new G4KineticTrackVector();
-     
+
      for (unsigned int aResult=0; aResult < result1->size(); aResult++)
      {
@@ -101,5 +105,5 @@
          theTotalResult->push_back(theNew);            
        }
-       else if(aTrack->Get4Momentum().t() - aTrack->Get4Momentum().mag()>80*MeV)
+       else if(aTrack->Get4Momentum().t() - aTrack->Get4Momentum().mag()>CaptureThreshold)
        {
 	 G4ReactionProduct * theNew = new G4ReactionProduct(aTrack->GetDefinition());
@@ -143,5 +147,5 @@
        theCurrentNucleon = theNucleus->GetNextNucleon();
      }   
-     
+ 
      if(!theDeExcitation)
      {
@@ -152,8 +156,9 @@
        G4double residualMass =  
                 G4ParticleTable::GetParticleTable()->GetIonTable()->GetIonMass(aZ ,anA);
-       residualMass += exEnergy;
+         residualMass += exEnergy;
+
        G4LorentzVector exciton4Momentum(exciton3Momentum, 
                                         std::sqrt(exciton3Momentum.mag2()+residualMass*residualMass));
-     
+    
        anInitialState.SetA(anA);
        anInitialState.SetZ(aZ);
@@ -167,5 +172,5 @@
        const G4Fragment aFragment(anInitialState);
        G4ReactionProductVector * aPreResult = theDeExcitation->DeExcite(aFragment);
-//       G4ReactionProductVector * aPreResult = new G4ReactionProductVector;
+
        // fill pre-compound part into the result, and return
        for(unsigned int ll=0; ll<aPreResult->size(); ll++)
@@ -186,2 +191,9 @@
    }
   
+G4double G4GeneratorPrecompoundInterface::SetCaptureThreshold(G4double value)
+{
+    G4double old=CaptureThreshold;
+    CaptureThreshold=value;
+    return old;
+
+}
Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4KM_NucleonEqRhs.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4KM_NucleonEqRhs.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4KM_NucleonEqRhs.cc	(revision 1196)
@@ -38,5 +38,4 @@
 
 #include "G4KM_NucleonEqRhs.hh"
-#include "G4NucleiPropertiesTable.hh"
 #include "G4VNuclearDensity.hh"
 
Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4NeutronField.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4NeutronField.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4NeutronField.cc	(revision 1196)
@@ -37,5 +37,4 @@
 // -------------------------------------------------------------------
 #include "G4NeutronField.hh"
-#include "G4NucleiPropertiesTable.hh"
 #include "G4VNuclearDensity.hh"
 #include "G4FermiMomentum.hh"
Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4ProtonField.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4ProtonField.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4ProtonField.cc	(revision 1196)
@@ -37,5 +37,4 @@
 // -------------------------------------------------------------------
 #include "G4ProtonField.hh"
-#include "G4NucleiPropertiesTable.hh"
 #include "G4VNuclearDensity.hh"
 #include "G4FermiMomentum.hh"
Index: trunk/source/processes/hadronic/models/binary_cascade/src/G4RKPropagation.cc
===================================================================
--- trunk/source/processes/hadronic/models/binary_cascade/src/G4RKPropagation.cc	(revision 962)
+++ trunk/source/processes/hadronic/models/binary_cascade/src/G4RKPropagation.cc	(revision 1196)
@@ -252,5 +252,5 @@
 {
 //  reset momentum transfer to field
-    theMomentumTranfer=0;
+    theMomentumTranfer=G4ThreeVector(0,0,0);
 
 // Loop over tracks
@@ -262,4 +262,5 @@
     G4KineticTrack * kt = *i;
     G4int encoding = kt->GetDefinition()->GetPDGEncoding();
+
     std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding);
 
@@ -327,7 +328,5 @@
 //     GetField = Barrier + FermiPotential
 	G4double newE = kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition());
-//      G4cout << " enter nucleus, E out/in: " << kt->GetTrackingMomentum().e() << " / " << newE <<G4endl;
-//      G4cout << " the Field "<< currentField->GetField(kt->GetPosition()) << " "<< kt->GetPosition()<<G4endl;
-//      G4cout << " the particle "<<kt->GetDefinition()->GetParticleName()<<G4endl;
+
 	if(newE <= kt->GetActualMass())  // the particle cannot enter the nucleus
 	{
@@ -350,12 +349,15 @@
 	kt->SetTrackingMomentum(new4Mom);
         kt->SetState(G4KineticTrack::inside);
-//     G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / "
-//    	   << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / "
-//	   << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition())
-//	   << G4endl
-//	   << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())"
-//	   << (*theFieldMap)[encoding]->GetBarrier() << " / "
-//	   << (*theFieldMap)[encoding]->GetField(0.9999*kt->GetPosition())
-//	   << G4endl;
+
+/*
+     G4cout <<" Enter Nucleus - E/Field/Sum: " <<kt->GetTrackingMomentum().e() << " / "
+    	   << (*theFieldMap)[encoding]->GetField(kt->GetPosition()) << " / "
+	   << kt->GetTrackingMomentum().e()-currentField->GetField(kt->GetPosition())
+	   << G4endl
+	   << " Barrier / field just inside nucleus (0.9999*kt->GetPosition())"
+	   << (*theFieldMap)[encoding]->GetBarrier() << " / "
+	   << (*theFieldMap)[encoding]->GetField(0.9999*kt->GetPosition())
+	   << G4endl;
+*/
        }
     }
@@ -396,5 +398,5 @@
  	<< G4endl << currentField->GetField(kt->GetPosition())<< G4endl
  	<< kt->GetTrackingMomentum()
- //	<< G4endl;
+ 	<< G4endl
 	<< "delta p " << momold-kt->GetTrackingMomentum() << G4endl
 	<< "del pos " << posold-kt->GetPosition()
@@ -483,6 +485,4 @@
       }
 
-
-
   }
 
@@ -502,5 +502,5 @@
 //----------------------------------------------------------------------------
 {
-    theMomentumTranfer=0;
+    theMomentumTranfer=G4ThreeVector(0,0,0);
 //    G4cout <<"Stepper input"<<kt->GetTrackingMomentum()<<G4endl;
 // create the integrator stepper
@@ -538,9 +538,9 @@
     }
 /*
- *      G4cout <<" E/Field/Sum be4 : " <<mom.e() << " / "
- *     	   << (*theFieldMap)[encoding]->GetField(pos) << " / "
- * 	   << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
- * 	   << G4endl;
- */
+       G4cout <<" E/Field/Sum be4 : " <<mom.e() << " / "
+      	   << (*theFieldMap)[encoding]->GetField(pos) << " / "
+  	   << mom.e()+(*theFieldMap)[encoding]->GetField(pos)
+  	   << G4endl;
+*/
 
 // Correct for momentum ( thus energy) transfered to nucleus, boost particle into moving nuclues frame.
