Changeset 1196 for trunk/source/processes/hadronic/models/binary_cascade/src

Timestamp:

Nov 25, 2009, 5:13:58 PM (15 years ago)

Author:

garnier

Message:

update CVS release candidate geant4.9.3.01

Location:

trunk/source/processes/hadronic/models/binary_cascade/src

Files:

• G4Absorber.cc (modified) (1 diff)
• G4BinaryCascade.cc (modified) (23 diffs)
• G4BinaryLightIonReaction.cc (modified) (3 diffs)
• G4GeneratorPrecompoundInterface.cc (modified) (7 diffs)
• G4KM_NucleonEqRhs.cc (modified) (1 diff)
• G4NeutronField.cc (modified) (1 diff)
• G4ProtonField.cc (modified) (1 diff)
• G4RKPropagation.cc (modified) (8 diffs)

trunk/source/processes/hadronic/models/binary_cascade/src/G4Absorber.cc

@@ -56 +56 @@
   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())
       {

trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryCascade.cc

@@ -97 +97 @@
   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 +107 @@
   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 +155 @@
   eventcounter++;
   if(getenv("BCDEBUG") ) G4cerr << " ######### Binary Cascade Reaction number starts ######### "<<eventcounter<<G4endl;
+
   G4LorentzVector initial4Momentum = aTrack.Get4Momentum();
   G4ParticleDefinition * definition = const_cast<G4ParticleDefinition *>(aTrack.GetDefinition());
@@ -238 +239 @@
   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 +249 @@
     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 +313 @@
   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 +360 @@
   G4bool haveProducts = false;
   G4int collisionCount=0;
-  theMomentumTransfer=0;
+  theMomentumTransfer=G4ThreeVector(0,0,0);
   while(theCollisionMgr->Entries() > 0)
   {
@@ -394 +391 @@
        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 +536 @@
   }
 
-
 // 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 +694 @@
          aNew->SetTotalEnergy(kt->Get4Momentum().e());
          Ekinout += aNew->GetKineticEnergy();
+         pSumBic += kt->Get4Momentum();
          if(kt->IsParticipant()) 
          {
@@ -655 +714 @@
              } else
              {
-                G4cout << "final un stable : ";
+                G4cout << "final non stable : ";
              }
              G4cout <<kt->GetDefinition()->GetParticleName()<< G4endl;
@@ -666 +725 @@
   }
 // add precompound products to products
+  G4LorentzVector pSumPreco(0), pPreco(0);
   if ( precompoundProducts )
   {
@@ -673 +733 @@
 // 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 +741 @@
          G4cout << " pProduct aft boost " <<pProduct << G4endl;
 #endif
+         pSumPreco += pProduct;
          (*j)->SetTotalEnergy(pProduct.e());
          (*j)->SetMomentum(pProduct.vect());
@@ -685 +747 @@
          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 +803 @@
      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 +994 @@
 //----------------------------------------------------------------------------
 {
+    G4RKPropagation * RKprop=(G4RKPropagation *)thePropagator;
     if ( secondary->GetTrackingMomentum().mag2() < -1.*eV )
     {
@@ -952 +1018 @@
     }
 
+    //  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 +1094 @@
   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 +1117 @@
 
   #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 +1207 @@
    }
 
-#ifdef debug_BIC_ApplyCollisionXX
-  DebugApplyCollsion(collision, products);
+#ifdef debug_BIC_ApplyCollision
+  DebugApplyCollision(collision, products);
 #endif
 
@@ -1231 +1315 @@
 }
 
-//----------------------------------------------------------------------------
-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 +2698 @@
 }
 
-
-
-
-
-
-
-
+//----------------------------------------------------------------------------
+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;
+
+}
+
+//----------------------------------------------------------------------------

trunk/source/processes/hadronic/models/binary_cascade/src/G4BinaryLightIonReaction.cc

@@ -38 +38 @@
 #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 +122 @@
 //      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 +406 @@
         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() )
       {

trunk/source/processes/hadronic/models/binary_cascade/src/G4GeneratorPrecompoundInterface.cc

@@ -31 +31 @@
 // HPW, 10DEC 98, the decay part originally written by Gunter Folger in his FTF-test-program.
 //
+   
+    G4GeneratorPrecompoundInterface::G4GeneratorPrecompoundInterface()
+    : CaptureThreshold(80*MeV)
+    {}
       
    G4HadFinalState* G4GeneratorPrecompoundInterface::
@@ -51 +55 @@
      result1=theSecondaries;
      result=new G4KineticTrackVector();
-     
+
      for (unsigned int aResult=0; aResult < result1->size(); aResult++)
      {
@@ -101 +105 @@
          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 +147 @@
        theCurrentNucleon = theNucleus->GetNextNucleon();
      }   
-     
+ 
      if(!theDeExcitation)
      {
@@ -152 +156 @@
        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 +172 @@
        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 +191 @@
    }
   
+G4double G4GeneratorPrecompoundInterface::SetCaptureThreshold(G4double value)
+{
+    G4double old=CaptureThreshold;
+    CaptureThreshold=value;
+    return old;
+
+}

trunk/source/processes/hadronic/models/binary_cascade/src/G4KM_NucleonEqRhs.cc

@@ -38 +38 @@
 
 #include "G4KM_NucleonEqRhs.hh"
-#include "G4NucleiPropertiesTable.hh"
 #include "G4VNuclearDensity.hh"
 

trunk/source/processes/hadronic/models/binary_cascade/src/G4NeutronField.cc

@@ -37 +37 @@
 // -------------------------------------------------------------------
 #include "G4NeutronField.hh"
-#include "G4NucleiPropertiesTable.hh"
 #include "G4VNuclearDensity.hh"
 #include "G4FermiMomentum.hh"

trunk/source/processes/hadronic/models/binary_cascade/src/G4ProtonField.cc

@@ -37 +37 @@
 // -------------------------------------------------------------------
 #include "G4ProtonField.hh"
-#include "G4NucleiPropertiesTable.hh"
 #include "G4VNuclearDensity.hh"
 #include "G4FermiMomentum.hh"

trunk/source/processes/hadronic/models/binary_cascade/src/G4RKPropagation.cc

@@ -252 +252 @@
 {
 //  reset momentum transfer to field
-    theMomentumTranfer=0;
+    theMomentumTranfer=G4ThreeVector(0,0,0);
 
 // Loop over tracks
@@ -262 +262 @@
     G4KineticTrack * kt = *i;
     G4int encoding = kt->GetDefinition()->GetPDGEncoding();
+
     std::map <G4int, G4VNuclearField*, std::less<G4int> >::iterator fieldIter= theFieldMap->find(encoding);
 
@@ -327 +328 @@
 //     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 +349 @@
         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 +398 @@
         << G4endl << currentField->GetField(kt->GetPosition())<< G4endl
         << kt->GetTrackingMomentum()
- //     << G4endl;
+        << G4endl
         << "delta p " << momold-kt->GetTrackingMomentum() << G4endl
         << "del pos " << posold-kt->GetPosition()
@@ -483 +485 @@
       }
 
-
-
   }
 
@@ -502 +502 @@
 //----------------------------------------------------------------------------
 {
-    theMomentumTranfer=0;
+    theMomentumTranfer=G4ThreeVector(0,0,0);
 //    G4cout <<"Stepper input"<<kt->GetTrackingMomentum()<<G4endl;
 // create the integrator stepper
@@ -538 +538 @@
     }
 /*
- *      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.