Changeset 962 for trunk/source/processes/hadronic/models/neutron_hp/src/G4NeutronHPDiscreteTwoBody.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/G4NeutronHPDiscreteTwoBody.cc (modified) (6 diffs)

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

@@ -28 +28 @@
 // A prototype of the low energy neutron transport model.
 //
+//080612 Bug fix contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2,3
+//080709 Bug fix Sampling Legendre expansion by T. Koi   
+//
 #include "G4NeutronHPDiscreteTwoBody.hh"
 #include "G4Gamma.hh"
@@ -92 +95 @@
      if(theCoeff[it].GetRepresentation()==0)
      {
+//TK Legendre expansion
        G4NeutronHPLegendreStore theStore(1);
        theStore.SetCoeff(0, theCoeff);
        theStore.SetManager(theManager);
-       cosTh = theStore.SampleMax(anEnergy);
+       //cosTh = theStore.SampleMax(anEnergy);
+       //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #3
+       cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
      }
      else if(theCoeff[it].GetRepresentation()==12) // means LINLIN
@@ -136 +142 @@
        if(theCoeff[it].GetRepresentation()==0)
        {
+//TK Legendre expansion
          G4NeutronHPLegendreStore theStore(2);
          theStore.SetCoeff(0, &(theCoeff[it-1]));
@@ -142 +149 @@
          aManager.Init(theManager.GetScheme(it), 2);
          theStore.SetManager(aManager);
-         cosTh = theStore.SampleMax(anEnergy);
+         //cosTh = theStore.SampleMax(anEnergy);
+//080709 TKDB
+         cosTh = theStore.SampleDiscreteTwoBody(anEnergy);
        }
        else if(theCoeff[it].GetRepresentation()==12) // LINLIN
@@ -198 +207 @@
          cosTh = theStore.Sample();
        }
-       else if(theCoeff[it].GetRepresentation()==14) 
+       else if(theCoeff[it].GetRepresentation()==14) //TK LOG_LIN
        {
          G4NeutronHPVector theBuff1;
@@ -266 +275 @@
 // now get the energy from kinematics and Q-value.
 
-   G4double restEnergy = anEnergy+GetQValue();
+   //G4double restEnergy = anEnergy+GetQValue();
    
 // assumed to be in CMS @@@@@@@@@@@@@@@@@
 
-   G4double residualMass =   GetTarget()->GetMass() + GetNeutron()->GetMass()
-                           - result->GetMass() - GetQValue();
-   G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
+   //080612TK contribution from Benoit Pirard and Laurent Desorgher (Univ. Bern) #2
+   //G4double residualMass =   GetTarget()->GetMass() + GetNeutron()->GetMass()
+   //                        - result->GetMass() - GetQValue();
+   //G4double kinE = restEnergy/(1+result->GetMass()/residualMass); // non relativistic @@
+   G4double A1     =  GetTarget()->GetMass()/GetNeutron()->GetMass(); 
+   G4double A1prim =  result->GetMass()/GetNeutron()->GetMass();
+   G4double E1     =  (A1+1)*(A1+1)/A1/A1*anEnergy; 
+   G4double kinE = (A1+1-A1prim)/(A1+1)/(A1+1)*(A1*E1+(1+A1)*GetQValue());
+
    result->SetKineticEnergy(kinE); // non relativistic @@
    G4double phi = twopi*G4UniformRand();