Changeset 962 for trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticCompFS.cc

Timestamp:

Apr 6, 2009, 12:30:29 PM (15 years ago)

Author:

garnier

Message:

update processes

File:

• trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPInelasticCompFS.cc (modified) (15 diffs)

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

@@ -30 +30 @@
 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
 // 070606 bug fix and migrate to enable to Partial cases by T. Koi 
+// 080603 bug fix for Hadron Hyper News #932 by T. Koi 
+// 080612 bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4,6
+// 080717 bug fix of calculation of residual momentum by T. Koi
+// 080801 protect negative avalable energy by T. Koi
+//        introduce theNDLDataA,Z which has A and Z of NDL data by T. Koi
+// 081024 G4NucleiPropertiesTable:: to G4NucleiProperties::
 //
 #include "G4NeutronHPInelasticCompFS.hh"
 #include "G4Nucleus.hh"
-#include "G4NucleiPropertiesTable.hh"
+#include "G4NucleiProperties.hh"
 #include "G4He3.hh"
 #include "G4Alpha.hh"
@@ -74 +80 @@
   theBaseA = aFile.GetA();
   theBaseZ = aFile.GetZ();
+   theNDLDataA = (int)aFile.GetA();
+   theNDLDataZ = aFile.GetZ();
   if(!dbool || ( Z<2.5 && ( std::abs(theBaseZ - Z)>0.0001 || std::abs(theBaseA - A)>0.0001)))
   {
@@ -185 +193 @@
 }
 
+
+                                                                                                       //n,p,d,t,he3,a
 void G4NeutronHPInelasticCompFS::CompositeApply(const G4HadProjectile & theTrack, G4ParticleDefinition * aDefinition)
 {
 
-    //G4cout << "G4NeutronHPInelasticCompFS::CompositeApply " << G4endl; 
 // prepare neutron
     theResult.Clear();
@@ -201 +210 @@
     for(i=0; i<50; i++)
     { if(theXsection[i] != 0) { break; } } 
+
     G4double targetMass=0;
     G4double eps = 0.0001;
-    targetMass = ( G4NucleiPropertiesTable::GetNuclearMass(static_cast<G4int>(theBaseZ+eps), static_cast<G4int>(theBaseA+eps))) /
+    targetMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(theBaseA+eps), static_cast<G4int>(theBaseZ+eps))) /
                    G4Neutron::Neutron()->GetPDGMass();
 //    if(theEnergyAngData[i]!=0)
@@ -220 +230 @@
     G4double residualZ = theBaseZ - aDefinition->GetPDGCharge();
     G4double residualA = theBaseA - aDefinition->GetBaryonNumber()+1;
-    residualMass = ( G4NucleiPropertiesTable::GetNuclearMass(static_cast<G4int>(residualZ+eps), static_cast<G4int>(residualA+eps)) ) /
+    residualMass = ( G4NucleiProperties::GetNuclearMass(static_cast<G4int>(residualA+eps), static_cast<G4int>(residualZ+eps)) ) /
                      G4Neutron::Neutron()->GetPDGMass();
 
@@ -230 +240 @@
   
 // select exit channel for composite FS class.
-    G4int it = SelectExitChannel(eKinetic);
+    G4int it = SelectExitChannel( eKinetic );
    
 // set target and neutron in the relevant exit channel
@@ -241 +251 @@
     G4double availableEnergy = theNeutron.GetKineticEnergy() + theNeutron.GetMass() - aHadron.GetMass() +
                              (targetMass - residualMass)*G4Neutron::Neutron()->GetPDGMass();
+//080730c
+    if ( availableEnergy < 0 )
+    {
+       //G4cout << "080730c Adjust availavleEnergy " << G4endl; 
+       availableEnergy = 0; 
+    }
     G4int nothingWasKnownOnHadron = 0;
     G4int dummy;
     G4double eGamm = 0;
     G4int iLevel=it-1;
-//  TK debug 070530  (without photon has it = 0)
-    //if(50==it) 
-    if( 0 == it ) 
-    {
+
+//  TK without photon has it = 0
+    if( 50 == it ) 
+    {
+
+//    TK Excitation level is not determined
       iLevel=-1;
       aHadron.SetKineticEnergy(availableEnergy*residualMass*G4Neutron::Neutron()->GetPDGMass()/
                                (aHadron.GetMass()+residualMass*G4Neutron::Neutron()->GetPDGMass()));
 
-      //aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*
-      //                  std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
-      //                            aHadron.GetMass()*aHadron.GetMass()));
-
+      aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*
+                        std::sqrt(aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-
+                                  aHadron.GetMass()*aHadron.GetMass()));
+
+/*
       G4double p2 = ( aHadron.GetTotalEnergy()*aHadron.GetTotalEnergy()-aHadron.GetMass()*aHadron.GetMass() );
       G4double p = 0.0;
@@ -264 +283 @@
       } 
       aHadron.SetMomentum(theNeutron.GetMomentum()*(1./theNeutron.GetTotalMomentum())*p );
+*/
 
     }
@@ -270 +290 @@
       while( iLevel!=-1 && theGammas.GetLevel(iLevel)==0 ) { iLevel--; }
     }
-    if(theAngularDistribution[it]!= 0)
-    {
-      if(theEnergyDistribution[it]!=0)
+
+
+    if ( theAngularDistribution[it] != 0 ) // MF4
+    {
+      if(theEnergyDistribution[it]!=0) // MF5
       {
         aHadron.SetKineticEnergy(theEnergyDistribution[it]->Sample(eKinetic, dummy));
@@ -304 +326 @@
       }
       theAngularDistribution[it]->SampleAndUpdate(aHadron);
-      if(theFinalStatePhotons[it] == 0)
+
+      if( theFinalStatePhotons[it] == 0 )
       {
 // TK comment Most n,n* eneter to this  
@@ -324 +347 @@
       }
     }
-    else if(theEnergyAngData[it] != 0)  
+    else if(theEnergyAngData[it] != 0) // MF6  
     {
       theParticles = theEnergyAngData[it]->Sample(eKinetic);
@@ -333 +356 @@
       nothingWasKnownOnHadron = 1;
     }
+
     //G4cout << "theFinalStatePhotons it " << it << G4endl;
     //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
-//  TK 070530
-    if ( it != 0 ) it = 50;  // it 50 has final state data for photon MF13 cross and MF14 ang
     //G4cout << "theFinalStatePhotons it " << it << G4endl;
     //G4cout << "theFinalStatePhotons[it] " << theFinalStatePhotons[it] << G4endl;
     //G4cout << "thePhotons " << thePhotons << G4endl;
-    if(theFinalStatePhotons[it]!=0) 
-    {
-      // the photon distributions are in the Nucleus rest frame.
+
+    if ( theFinalStatePhotons[it] != 0 ) 
+    {
+       // the photon distributions are in the Nucleus rest frame.
+       // TK residual rest frame
       G4ReactionProduct boosted;
       boosted.Lorentz(theNeutron, theTarget);
@@ -493 +517 @@
                                   ->GetIon(static_cast<G4int>(residualZ), static_cast<G4int>(residualA), 0));  
         theResidual.SetKineticEnergy(aHadron.GetKineticEnergy()*aHadron.GetMass()/theResidual.GetMass());
-        theResidual.SetMomentum(-1.*aHadron.GetMomentum());
+
+        //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #6
+        //theResidual.SetMomentum(-1.*aHadron.GetMomentum());
+        G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
+        theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
+
         theResidual.Lorentz(theResidual, -1.*theTarget);
         G4ThreeVector totalPhotonMomentum(0,0,0);
@@ -529 +558 @@
 //        cout << "Kinetic energy of the residual = "<<resiualKineticEnergy<<endl;
         theResidual.SetKineticEnergy(resiualKineticEnergy);
-        theResidual.SetMomentum(-1.*totalMomentum);
+
+        //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #4
+        //theResidual.SetMomentum(-1.*totalMomentum);
+        //G4ThreeVector incidentNeutronMomentum = theNeutron.GetMomentum();
+        //theResidual.SetMomentum(incidentNeutronMomentum - aHadron.GetMomentum());
+//080717 TK Comment still do NOT include photon's mometum which produce by thePhotons
+        theResidual.SetMomentum( theNeutron.GetMomentum() + theTarget.GetMomentum() - totalMomentum );
+
         theSec = new G4DynamicParticle;   
         theSec->SetDefinition(theResidual.GetDefinition());
@@ -549 +585 @@
       delete thePhotons;
     }
+
+//080721 
+   G4ParticleDefinition* targ_pd = G4ParticleTable::GetParticleTable()->GetIon ( (G4int)theBaseZ , (G4int)theBaseA , 0.0 );
+   G4LorentzVector targ_4p_lab ( theTarget.GetMomentum() , std::sqrt( targ_pd->GetPDGMass()*targ_pd->GetPDGMass() + theTarget.GetMomentum().mag2() ) );
+   G4LorentzVector proj_4p_lab = theTrack.Get4Momentum();
+   G4LorentzVector init_4p_lab = proj_4p_lab + targ_4p_lab;
+   adjust_final_state ( init_4p_lab ); 
+
 // clean up the primary neutron
     theResult.SetStatusChange(stopAndKill);