Changeset 1347 for trunk/source/processes/hadronic/models/neutron_hp/src

Timestamp:

Dec 22, 2010, 3:52:27 PM (14 years ago)

Author:

garnier

Message:

geant4 tag 9.4

Location:

trunk/source/processes/hadronic/models/neutron_hp/src

Files:

• G4NeutronHPCaptureFS.cc (modified) (5 diffs)
• G4NeutronHPDiscreteTwoBody.cc (modified) (7 diffs)
• G4NeutronHPFinalState.cc (modified) (5 diffs)
• G4NeutronHPInelasticBaseFS.cc (modified) (4 diffs)
• G4NeutronHPInelasticCompFS.cc (modified) (2 diffs)
• G4NeutronHPPhotonDist.cc (modified) (2 diffs)

trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPCaptureFS.cc

@@ -31 +31 @@
 // 26-January-07 Add G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION flag
 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
+// 101203 Bugzilla/Geant4 Problem 1155 Lack of residual in some case 
 //
 #include "G4NeutronHPCaptureFS.hh"
@@ -43 +44 @@
   G4HadFinalState * G4NeutronHPCaptureFS::ApplyYourself(const G4HadProjectile & theTrack)
   {
+
     G4int i;
     theResult.Clear();
@@ -108 +110 @@
     }
 
+
+
 // add them to the final state
 
@@ -114 +118 @@
     G4int nParticles = nPhotons;
     if(1==nPhotons) nParticles = 2;
+
+
+//Make at least one photon  
+//101203 TK
+    if ( nPhotons == 0 )
+    {
+       G4ReactionProduct * theOne = new G4ReactionProduct;
+       theOne->SetDefinition( G4Gamma::Gamma() ); 
+       G4double theta = pi*G4UniformRand();
+       G4double phi = twopi*G4UniformRand();
+       G4double sinth = std::sin(theta);
+       G4ThreeVector direction( sinth*std::cos(phi), sinth*std::sin(phi), std::cos(theta) );
+       theOne->SetMomentum( direction ) ;
+       thePhotons->push_back(theOne);
+       nPhotons++; // 0 -> 1
+    }
+//One photon case: energy set to Q-value 
+//101203 TK
+    if ( nPhotons == 1 )
+    {
+       G4ThreeVector direction = thePhotons->operator[](0)->GetMomentum().unit();
+       G4double Q = G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA), 0, static_cast<G4int>(theBaseZ))->GetPDGMass() + G4Neutron::Neutron()->GetPDGMass()
+         - G4ParticleTable::GetParticleTable()->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ))->GetPDGMass();
+       thePhotons->operator[](0)->SetMomentum( Q*direction );
+    } 
+//
 
     // back to lab system
@@ -156 +186 @@
     }
     delete thePhotons; 
+
+//101203TK
+    G4bool residual = false;
+    G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
+                                   ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
+    for ( G4int i = 0 ; i != theResult.GetNumberOfSecondaries() ; i++ )
+    {
+       if ( theResult.GetSecondary(i)->GetParticle()->GetDefinition() == aRecoil ) residual = true;
+    }
+
+    if ( residual == false )
+    {
+       G4ParticleDefinition * aRecoil = G4ParticleTable::GetParticleTable()
+                                        ->FindIon(static_cast<G4int>(theBaseZ), static_cast<G4int>(theBaseA+1), 0, static_cast<G4int>(theBaseZ));
+       G4int nNonZero = 0;
+       G4LorentzVector p_photons(0,0,0,0);
+       for ( G4int i = 0 ; i != theResult.GetNumberOfSecondaries() ; i++ )
+       {
+          p_photons += theResult.GetSecondary(i)->GetParticle()->Get4Momentum();
+          // To many 0 momentum photons -> Check PhotonDist 
+          if ( theResult.GetSecondary(i)->GetParticle()->Get4Momentum() > 0 ) nNonZero++;
+       }
+
+       // Can we include kinetic energy here?
+       G4double deltaE = ( theTrack.Get4Momentum().e() + theTarget.GetTotalEnergy() )
+                       - ( p_photons.e() + aRecoil->GetPDGMass() );
+
+//Add photons
+       if ( nPhotons - nNonZero > 0 ) 
+       {
+              //G4cout << "TKDB G4NeutronHPCaptureFS::ApplyYourself we will create additional " << nPhotons - nNonZero << " photons" << G4endl;
+          std::vector<G4double> vRand;
+          vRand.push_back( 0.0 );
+          for ( G4int i = 0 ; i != nPhotons - nNonZero - 1 ; i++ )
+          { 
+             vRand.push_back( G4UniformRand() );
+          }
+          vRand.push_back( 1.0 );
+          std::sort( vRand.begin(), vRand.end() );
+
+          std::vector<G4double> vEPhoton;
+          for ( G4int i = 0 ; i < (G4int)vRand.size() - 1 ; i++ )
+          {
+             vEPhoton.push_back( deltaE * ( vRand[i+1] - vRand[i] ) );
+          }
+          std::sort( vEPhoton.begin(), vEPhoton.end() );
+
+          for ( G4int i = 0 ; i < nPhotons - nNonZero - 1 ; i++ )
+          {
+             //Isotopic in LAB OK?
+             G4double theta = pi*G4UniformRand();
+             G4double phi = twopi*G4UniformRand();
+             G4double sinth = std::sin(theta);
+             G4double en = vEPhoton[i];
+             G4ThreeVector tempVector(en*sinth*std::cos(phi), en*sinth*std::sin(phi), en*std::cos(theta) );
+              
+             p_photons += G4LorentzVector ( tempVector, tempVector.mag() );
+             G4DynamicParticle * theOne = new G4DynamicParticle;
+             theOne->SetDefinition( G4Gamma::Gamma() );
+             theOne->SetMomentum( tempVector );
+             theResult.AddSecondary(theOne);
+          }
+
+//        Add last photon 
+          G4DynamicParticle * theOne = new G4DynamicParticle;
+          theOne->SetDefinition( G4Gamma::Gamma() );
+//        For better momentum conservation 
+          G4ThreeVector lastPhoton = -p_photons.vect().unit()*vEPhoton.back();
+          p_photons += G4LorentzVector( lastPhoton , lastPhoton.mag() );
+          theOne->SetMomentum( lastPhoton );
+          theResult.AddSecondary(theOne);
+       }
+
+//Add residual 
+       G4DynamicParticle * theOne = new G4DynamicParticle;
+       G4ThreeVector aMomentum = theTrack.Get4Momentum().vect() + theTarget.GetMomentum()
+                               - p_photons.vect();
+       theOne->SetDefinition(aRecoil);
+       theOne->SetMomentum( aMomentum );
+       theResult.AddSecondary(theOne);
+
+    }
+//101203TK END
+
 // clean up the primary neutron
     theResult.SetStatusChange(stopAndKill);

trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPDiscreteTwoBody.cc

@@ -30 +30 @@
 //080612 Bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2,3
 //080709 Bug fix Sampling Legendre expansion by T. Koi   
+//101110 Bug fix in MF=6, LAW=2 case; contribution from E. Mendoza, D. Cano-Ott (CIEMAT)
 //
 #include "G4NeutronHPDiscreteTwoBody.hh"
@@ -111 +112 @@
        for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
        {
-         theStore.SetX(i, theCoeff[it].GetCoeff(i));
-         theStore.SetY(i, theCoeff[it].GetCoeff(i));
+         //101110
+         //theStore.SetX(i, theCoeff[it].GetCoeff(i));
+         //theStore.SetY(i, theCoeff[it].GetCoeff(i));
+         theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
+         theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
          i++;
        }
@@ -125 +129 @@
        for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
        {
-         theStore.SetX(i, theCoeff[it].GetCoeff(i));
-         theStore.SetY(i, theCoeff[it].GetCoeff(i));
+         //101110
+         //theStore.SetX(i, theCoeff[it].GetCoeff(i));
+         //theStore.SetY(i, theCoeff[it].GetCoeff(i));
+         theStore.SetX(i/2, theCoeff[it].GetCoeff(i));
+         theStore.SetY(i/2, theCoeff[it].GetCoeff(i+1));
          i++;
        }
@@ -161 +168 @@
          for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
          {
-           theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
-           theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
+           //101110
+           //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
+           //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
+           theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
+           theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
            i++;
          }
@@ -171 +181 @@
          for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
          {
+           //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
+           //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
            theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
-           theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
+           theBuff2.SetY(i, theCoeff[it].GetCoeff(i+1));
            i++;
          }
@@ -215 +227 @@
          for(i=0;i<theCoeff[it-1].GetNumberOfPoly(); i++)
          {
-           theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
-           theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
+           //101110
+           //theBuff1.SetX(i, theCoeff[it-1].GetCoeff(i));
+           //theBuff1.SetY(i, theCoeff[it-1].GetCoeff(i));
+           theBuff1.SetX(i/2, theCoeff[it-1].GetCoeff(i));
+           theBuff1.SetY(i/2, theCoeff[it-1].GetCoeff(i+1));
            i++;
          }
@@ -226 +241 @@
          for(i=0;i<theCoeff[it].GetNumberOfPoly(); i++)
          {
-           theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
-           theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
+           //101110
+           //theBuff2.SetX(i, theCoeff[it].GetCoeff(i));
+           //theBuff2.SetY(i, theCoeff[it].GetCoeff(i));
+           theBuff2.SetX(i/2, theCoeff[it].GetCoeff(i));
+           theBuff2.SetY(i/2, theCoeff[it].GetCoeff(i+1));
            i++;
          }

trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPFinalState.cc

@@ -27 +27 @@
 //080721 Create adjust_final_state method by T. Koi
 //080801 Residual reconstruction with theNDLDataA,Z (A, Z, and momentum are adjusted) by T. Koi
+//101110 Set lower limit for gamma energy(1keV) by T. Koi
 
 #include "G4NeutronHPFinalState.hh"
@@ -37 +38 @@
 void G4NeutronHPFinalState::adjust_final_state ( G4LorentzVector init_4p_lab )
 {
+
+   G4double minimum_energy = 1*keV;
 
    if ( adjustResult != true ) return;
@@ -175 +178 @@
 
       // Adjust p
-      if ( dif_4p.v().mag() < 1*MeV )
+      //if ( dif_4p.v().mag() < 1*MeV )
+      if ( minimum_energy < dif_4p.v().mag() && dif_4p.v().mag() < 1*MeV )
       {
 
@@ -184 +188 @@
       else
       {
-         //G4cout << "HP_DB Difference in dif_p is too large (>1MeV) to adjust, so that give up tuning" << G4endl;
+         //G4cout << "HP_DB Difference in dif_p is too large (>1MeV) or too small(<1keV) to adjust, so that give up tuning" << G4endl;
       }
 
@@ -213 +217 @@
       nSecondaries += 2;
       G4double e1 = ( dif_4p.e() -dif_4p.v().mag() ) / 2;
-      G4double costh = 2.*G4UniformRand()-1.;
-      G4double phi = twopi*G4UniformRand();
-      G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi), 
-                         std::sin(std::acos(costh))*std::sin(phi),
-                         costh);
-      theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );    
-      theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );    
+      
+      if ( minimum_energy < e1 )  
+      {
+         G4double costh = 2.*G4UniformRand()-1.;
+         G4double phi = twopi*G4UniformRand();
+         G4ThreeVector dir( std::sin(std::acos(costh))*std::cos(phi), 
+                            std::sin(std::acos(costh))*std::sin(phi),
+                            costh);
+         theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , e1*dir ) );    
+         theResult.AddSecondary ( new G4DynamicParticle ( G4Gamma::Gamma() , -e1*dir ) );    
+      }
+      else
+      {
+         //G4cout << "HP_DB Difference is too small(<1keV) to adjust, so that neglect it" << G4endl;
+      }
 
    }

trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticBaseFS.cc

@@ -31 +31 @@
 //        Introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
 // 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
+// 101111 Add Special treatment for Be9(n,2n)Be8(2a) case by T. Koi
 //
 #include "G4NeutronHPInelasticBaseFS.hh"
@@ -360 +361 @@
   else if(theEnergyAngData!=0)
   {
+
     G4double theGammaEnergy = theEnergyAngData->GetTotalMeanEnergy();
     G4double anEnergy = boosted.GetKineticEnergy();
@@ -371 +373 @@
     G4double eBindHe3 = G4NucleiProperties::GetBindingEnergy(3,2);
     G4double eBindA = G4NucleiProperties::GetBindingEnergy(4,2);
+    G4int ia=0;
     for(i=0; i<tmpHadrons->size(); i++)
     {
@@ -396 +399 @@
       {
         eBindProducts+=eBindA;
-      }
-    }
+        ia++; 
+      }
+    }
+
     theGammaEnergy += eBindProducts;
+
+//101111 
+//Special treatment for Be9 + n -> 2n + Be8 -> 2n + a + a
+if ( (G4int)(theBaseZ+eps) == 4 && (G4int)(theBaseA+eps) == 9 )
+{
+   // This only valid for G4NDL3.13,,,
+   if ( std::abs( theNuclearMassDifference -   
+        ( G4NucleiProperties::GetBindingEnergy( 8 , 4 ) -
+        G4NucleiProperties::GetBindingEnergy( 9 , 4 ) ) ) < 1*keV 
+      && ia == 2 )
+   {
+      theGammaEnergy -= (2*eBindA);
+   }
+}
     
     G4ReactionProductVector * theOtherPhotons = 0;

trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticCompFS.cc

@@ -41 +41 @@
 //        add two_body_reaction
 // 100909 add safty 
+// 101111 add safty for _nat_ data case in Binary reaction, but break conservation  
 //
 #include "G4NeutronHPInelasticCompFS.hh"
@@ -681 +682 @@
    G4double AA = hadron->GetDefinition()->GetPDGMass() / proj->GetDefinition()->GetPDGMass(); 
    G4double E1 = proj->GetKineticEnergy();
-   G4double beta = std::sqrt ( A*(A+1-AA)/AA*(1+(1+A)/A*Q/E1) );
+
+// 101111
+// In _nat_ data (Q+E1) could become negative value, following line is safty for this case.
+   //if ( (Q+E1) < 0 ) 
+   if ( ( 1 + (1+A)/A*Q/E1 ) < 0 ) 
+   {
+// 1.0e-6 eV is additional safty for numeric precision 
+      Q = -( A/(1+A)*E1 ) + 1.0e-6*eV;
+   }
+
+   G4double beta = std::sqrt ( A*(A+1-AA)/AA*( 1 + (1+A)/A*Q/E1 ) );
    G4double gamma = AA/(A+1-AA)*beta;
    G4double E3 = AA/std::pow((1+A),2)*(beta*beta+1+2*beta*mu)*E1;

trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPPhotonDist.cc

@@ -39 +39 @@
 //        Contribution from Chao Zhang (Chao.Zhang@usd.edu) and Dongming Mei(Dongming.Mei@usd.edu)
 //        But it looks like never cause real effect in G4NDL3.13 (at least Natural elements) TK
+// 101111 Change warning message for "repFlag == 2 && isoFlag != 1" case 
 //
 // there is a lot of unused (and undebugged) code in this file. Kept for the moment just in case. @@
@@ -124 +125 @@
   if (isoFlag != 1)
   {
-if ( repFlag == 2 ) G4cout << "TKDB repFlag == 2 && isoFlag !=1  " << G4endl;
+if ( repFlag == 2 ) G4cout << "G4NeutronHPPhotonDist: repFlag == 2 && isoFlag != 1 is unexpected! If you use G4ND3.x, then please report to Geant4 Hyper News. Thanks." << G4endl;
     aDataFile >> tabulationType >> nDiscrete2 >> nIso;
 //080731