G4NeutronHPKallbachMannSyst Class Reference

#include <G4NeutronHPKallbachMannSyst.hh>


Public Member Functions

 G4NeutronHPKallbachMannSyst (G4double aCompoundFraction, G4double anIncidentEnergy, G4double anIncidentMass, G4double aProductEnergy, G4double aProductMass, G4double aResidualMass, G4int aResidualA, G4int aResidualZ, G4double aTargetMass, G4int aTargetA, G4int aTargetZ)
 ~G4NeutronHPKallbachMannSyst ()
G4double Sample (G4double anEnergy)
G4double Kallbach (G4double cosTh, G4double anEnergy)
G4double GetKallbachZero (G4double anEnergy)
G4double A (G4double anEnergy)
G4double SeparationEnergy (G4int Ac, G4int Nc, G4int AA, G4int ZA)

Detailed Description

Definition at line 34 of file G4NeutronHPKallbachMannSyst.hh.

Constructor & Destructor Documentation

G4NeutronHPKallbachMannSyst::G4NeutronHPKallbachMannSyst ( G4double  aCompoundFraction,
G4double  anIncidentEnergy,
G4double  anIncidentMass,
G4double  aProductEnergy,
G4double  aProductMass,
G4double  aResidualMass,
G4int  aResidualA,
G4int  aResidualZ,
G4double  aTargetMass,
G4int  aTargetA,
G4int  aTargetZ 
) [inline]

Definition at line 38 of file G4NeutronHPKallbachMannSyst.hh.

00043   {
00044     theCompoundFraction = aCompoundFraction;
00045     theIncidentEnergy = anIncidentEnergy;
00046     theIncidentMass = anIncidentMass;
00047     theProductEnergy = aProductEnergy;
00048     theProductMass = aProductMass;
00049     theResidualMass = aResidualMass;
00050     theResidualA = aResidualA;
00051     theResidualZ = aResidualZ;
00052     theTargetMass = aTargetMass;
00053     theTargetA = aTargetA;
00054     theTargetZ = aTargetZ;
00055   }

G4NeutronHPKallbachMannSyst::~G4NeutronHPKallbachMannSyst (  )  [inline]

Definition at line 57 of file G4NeutronHPKallbachMannSyst.hh.

00057 {};

Member Function Documentation

G4double G4NeutronHPKallbachMannSyst::A ( G4double  anEnergy  ) 

Definition at line 83 of file G4NeutronHPKallbachMannSyst.cc.

References C1, C3, and SeparationEnergy().

Referenced by GetKallbachZero(), and Kallbach().

00084 { 
00085   G4double result;
00086   G4double C1 = 0.04/MeV;
00087   G4double C2 = 1.8E-6/(MeV*MeV*MeV);
00088   G4double C3 = 6.7E-7/(MeV*MeV*MeV*MeV);
00089   
00090     G4double epsa = anEnergy*theTargetMass/(theTargetMass+theIncidentMass);
00091     G4int Ac = theTargetA+1; 
00092     G4int Nc = Ac - theTargetZ;
00093     G4int AA = theTargetA;
00094     G4int ZA = theTargetZ;
00095     G4double ea = epsa+SeparationEnergy(Ac, Nc, AA, ZA);
00096     G4double Et1 = 130*MeV;
00097     G4double R1 = std::min(ea, Et1);
00098     // theProductEnergy is still in CMS!!!
00099     G4double epsb = theProductEnergy*(theProductMass+theResidualMass)/theResidualMass;
00100     G4int AB = theResidualA;
00101     G4int ZB = theResidualZ;
00102     G4double eb = epsb+SeparationEnergy(Ac, Nc, AB, ZB );
00103   G4double X1 = R1*eb/ea; 
00104     G4double Et3 = 41*MeV;
00105     G4double R3 = std::min(ea, Et3);
00106   G4double X3 = R3*eb/ea;
00107   G4double Ma = 1;
00108   G4double mb(0);
00109   G4int productA = theTargetA+1-theResidualA;
00110   G4int productZ = theTargetZ-theResidualZ;
00111   if(productZ==0)
00112   {
00113     mb = 0.5;
00114   }
00115   else if(productZ==1)
00116   {
00117     mb = 1;
00118   }
00119   else if(productZ==2)
00120   {
00121     mb = 2;
00122     if(productA==3) mb=1;
00123   }
00124   else
00125   {
00126     throw G4HadronicException(__FILE__, __LINE__, "Severe error in the sampling of Kallbach-Mann Systematics");
00127   }
00128   
00129   result = C1*X1 + C2*std::pow(X1, 3.) + C3*Ma*mb*std::pow(X3, 4.);
00130   return result;
00131 }

G4double G4NeutronHPKallbachMannSyst::GetKallbachZero ( G4double  anEnergy  ) 

Definition at line 71 of file G4NeutronHPKallbachMannSyst.cc.

References A().

Referenced by Sample().

00072 {
00073   G4double result;
00074   if ( theCompoundFraction == 1 ) 
00075   { 
00076      //G4cout << "080730b Adjust theCompoundFraction " << G4endl;
00077      theCompoundFraction *= (1-1.0e-15);   
00078   } 
00079   result = 0.5 * (1./A(anEnergy)) * std::log((1-theCompoundFraction)/(1+theCompoundFraction));
00080   return result;
00081 }

G4double G4NeutronHPKallbachMannSyst::Kallbach ( G4double  cosTh,
G4double  anEnergy 
)

Definition at line 61 of file G4NeutronHPKallbachMannSyst.cc.

References A().

Referenced by Sample().

00062 {
00063   // Kallbach-Mann systematics without normalization.
00064   G4double result;
00065   G4double theX = A(anEnergy)*cosTh;
00066   result = 0.5*(std::exp( theX)*(1+theCompoundFraction)
00067                +std::exp(-theX)*(1-theCompoundFraction));
00068   return result;
00069 }

G4double G4NeutronHPKallbachMannSyst::Sample ( G4double  anEnergy  ) 

Definition at line 37 of file G4NeutronHPKallbachMannSyst.cc.

References G4UniformRand, GetKallbachZero(), and Kallbach().

Referenced by G4NeutronHPContAngularPar::Sample().

00038 {
00039   G4double result;
00040   
00041   G4double zero = GetKallbachZero(anEnergy);
00042   if(zero>1) zero=1.;
00043   if(zero<-1)zero=-1.;
00044   G4double max = Kallbach(zero, anEnergy);
00045   double upper = Kallbach(1., anEnergy);
00046   double lower = Kallbach(-1., anEnergy);
00047   if(upper>max) max=upper;
00048   if(lower>max) max=lower;
00049   G4double value, random;
00050   do
00051   {
00052     result = 2.*G4UniformRand()-1;  
00053     value = Kallbach(result, anEnergy)/max;
00054     random = G4UniformRand();
00055   }
00056   while(random>value);
00057   
00058   return result;
00059 }

G4double G4NeutronHPKallbachMannSyst::SeparationEnergy ( G4int  Ac,
G4int  Nc,
G4int  AA,
G4int  ZA 
)

Definition at line 133 of file G4NeutronHPKallbachMannSyst.cc.

Referenced by A().

00134 {
00135   G4double result;
00136   G4int NA = AA-ZA;
00137   G4int Zc = Ac-Nc;
00138   result = 15.68*(Ac-AA);
00139   result += -28.07*((Nc-Zc)*(Nc-Zc)/Ac - (NA-ZA)*(NA-ZA)/AA);
00140   result += -18.56*(std::pow(G4double(Ac), 2./3.) - std::pow(G4double(AA), 2./3.));
00141   result +=  33.22*((Nc-Zc)*(Nc-Zc)/std::pow(G4double(Ac), 4./3.) - (NA-ZA)*(NA-ZA)/std::pow(G4double(AA), 4./3.));
00142   result += -0.717*(Zc*Zc/std::pow(G4double(Ac),1./3.)-ZA*ZA/std::pow(G4double(AA),1./3.));
00143   result +=  1.211*(Zc*Zc/Ac-ZA*ZA/AA);
00144   G4double totalBinding(0);
00145   G4int productA = theTargetA+1-theResidualA;
00146   G4int productZ = theTargetZ-theResidualZ;
00147   if(productZ==0&&productA==1) totalBinding=0;
00148   if(productZ==1&&productA==1) totalBinding=0;
00149   if(productZ==1&&productA==2) totalBinding=2.22;
00150   if(productZ==1&&productA==3) totalBinding=8.48;
00151   if(productZ==2&&productA==3) totalBinding=7.72;
00152   if(productZ==2&&productA==4) totalBinding=28.3;
00153   result += -totalBinding;
00154   result *= MeV;
00155   return result;
00156 }

The documentation for this class was generated from the following files:
Generated on Mon May 27 17:52:39 2013 for Geant4 by  doxygen 1.4.7