Changeset 1228 for trunk/source/processes/hadronic/models/binary_cascade

Timestamp:

Jan 8, 2010, 11:56:51 AM (14 years ago)

Author:

garnier

Message:

update geant4.9.3 tag

Location:

trunk/source/processes/hadronic/models/binary_cascade

Files:

• History (modified) (1 diff)
• include/G4BinaryCascade.hh (modified) (1 diff)
• src/G4BinaryCascade.cc (modified) (10 diffs)

trunk/source/processes/hadronic/models/binary_cascade/History

@@ -13 +13 @@
      * Please list in reverse chronological order (last date on top)
      ---------------------------------------------------------------
+
+4-Dec-2009, Gunter Folger                had-binary-V09-02-06
+- Bug fix in G4BinaryCascade::ApplyCollision; decay products outside
+   nucleus were nevertheless counted in currentZ/A as if these were within
+   nucleus.
 
 13-Nov-2009, Gunter Folger               had-binary-V09-02-05

trunk/source/processes/hadronic/models/binary_cascade/include/G4BinaryCascade.hh

@@ -1 @@
-
+//
 // ********************************************************************
 // * License and Disclaimer                                           *

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

@@ -787 +787 @@
 // 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 +852 @@
       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 +1081 @@
   G4int initialBaryon(0);
   G4int initialCharge(0);
-  G4double initial_Efermi(0);
   G4LorentzVector mom4Primary(0);
   
@@ -1090 +1089 @@
   }
 
-// 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 +1142 @@
      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 +1174 @@
   }
 
-   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 +1222 @@
     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 +1253 @@
        }
        
-//       G4cout << " PDGcode, state " << (*i)->GetDefinition()->GetPDGEncoding() << " " << (*i)->GetState()<<G4endl;
+       //G4cout << " PDGcode, state " << (*i)->GetDefinition()->GetPDGEncoding() << " " << (*i)->GetState()<<G4endl;
 
     }   
@@ -1269 +1275 @@
   currentA += finalBaryon-initialBaryon;
   currentZ += finalCharge-initialCharge;
-        //G4cout << " currentA, Z aft: " << currentA << " " << currentZ << G4endl;
+        //G4cout << " ApplyCollision currentA, Z aft: " << currentA << " " << currentZ << G4endl;
   
   G4KineticTrackVector oldSecondaries;
@@ -2351 +2357 @@
      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)
      {