GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
INukeOsetFormula.cxx
Go to the documentation of this file.
1 #include "INukeOsetFormula.h"
2 
3 // constants
4 const double INukeOsetFormula :: fCouplingConstant = 0.36 * 4.0 * kPi;
5 
6 const double INukeOsetFormula :: fNormalDensity = 0.17; // fm-3
7 
8 const double INukeOsetFormula :: fNormFactor = 197.327 * 197.327 * 10.0;
9 
10 // particles mass
11 
12 const double INukeOsetFormula :: fNucleonMass = kNucleonMass * 1000.0; // [MeV]
13 const double INukeOsetFormula :: fNucleonMass2 = fNucleonMass * fNucleonMass;
14 
15 const double INukeOsetFormula :: fDeltaMass = 1232.0; // MeV
16 
17 // s-wave parametrization (see sec. 3.1)
18 
19 const double INukeOsetFormula :: fCoefSigma[3] = {-0.01334, 0.06889, 0.19753};
20 const double INukeOsetFormula :: fCoefB[3] = {-0.01866, 0.06602, 0.21972};
21 const double INukeOsetFormula :: fCoefD[3] = {-0.08229, 0.37062,-0.03130};
22 const double INukeOsetFormula :: ImB0 = 0.035;
23 
24 //! delta parametrization coefficients (NuclPhys A468 (1987) 631-652)
25 const double INukeOsetFormula :: fCoefCQ[3] = { -5.19, 15.35, 2.06};
26 //! delta parametrization coefficients (NuclPhys A468 (1987) 631-652)
27 const double INukeOsetFormula :: fCoefCA2[3] = { 1.06, -6.64, 22.66};
28 //! delta parametrization coefficients (NuclPhys A468 (1987) 631-652)
29 const double INukeOsetFormula :: fCoefCA3[3] = {-13.46, 46.17, -20.34};
30 //! delta parametrization coefficients (NuclPhys A468 (1987) 631-652)
31 const double INukeOsetFormula :: fCoefAlpha[3] = {0.382, -1.322, 1.466};
32 //! delta parametrization coefficients (NuclPhys A468 (1987) 631-652)
33 const double INukeOsetFormula :: fCoefBeta[3] = {-0.038, 0.204, 0.613};
34 
35 using namespace osetUtils;
36 
37 //! @f$ p_{F} = \left(\frac{3}{2}\pi\rho\right)^{1/3} @f$
38 void INukeOsetFormula :: setNucleus (const double &density)
39 {
40  static const double constFactor = 3.0 / 2.0 * kPi * kPi;
41 
42  fNuclearDensity = density;
43  fFermiMomentum = pow (constFactor * fNuclearDensity, 1.0 / 3.0) * 197.327;
44  fFermiEnergy = sqrt (fFermiMomentum*fFermiMomentum + fNucleonMass2);
45 }
46 
47 /*! <ul>
48  * <li> set pointer to proper mass (charged vs neutral)
49  * <li> calculate energy, momentum (in LAB and CMS) and invariant mass
50  * </ul>
51  */
52 void INukeOsetFormula :: setKinematics (const double &pionTk, const bool &isPi0)
53 {
54  if (isPi0)
55  {
56  fPionMass = PDGLibrary::Instance()->Find(kPdgPi0)->Mass() * 1000.0; // [MeV]
57  fPionMass2 = fPionMass * fPionMass;
58  }
59  else
60  {
61  fPionMass = PDGLibrary::Instance()->Find(kPdgPiP)->Mass() * 1000.0; // [MeV]
62  fPionMass2 = fPionMass * fPionMass;
63  }
64 
65  fPionKineticEnergy = pionTk;
66  fPionEnergy = fPionKineticEnergy + fPionMass;
67  fPionMomentum = sqrt (fPionEnergy * fPionEnergy - fPionMass2);
68 
69  // assuming nucleon at rest
70  fInvariantMass = sqrt (fPionMass2 + 2.0 * fNucleonMass * fPionEnergy +
71  fNucleonMass2);
72 
73  fMomentumCMS = fPionMomentum * fNucleonMass / fInvariantMass;
74  fMomentumCMS2 = fMomentumCMS * fMomentumCMS;
75 }
76 
77 /*! <ul>
78  * <li> calculale Delta half width (in nuclear matter)
79  * <li> calculalte Delta self energy
80  * <li> calculalte Delta propagator
81  * </ul>
82  */
84 {
85  static const double constFactor = 1.0 / 12.0 / kPi;
86 
87  fCouplingFactor = fCouplingConstant / fPionMass2; // (f*/m_pi)^2, e.g. eq. 2.6
88 
89  // see eq. 2.7 and sec. 2.3
90  fReducedHalfWidth = constFactor * fCouplingFactor * fNucleonMass * fMomentumCMS *
91  fMomentumCMS2 / fInvariantMass * deltaReduction ();
92 
93  setSelfEnergy (); // calculate qel and absorption part of delta self-energy
94 
95  // real and imaginary part of delta denominator (see eq. 2.23)
96  const double ReDelta = fInvariantMass - fDeltaMass;
97  const double ImDelta = fReducedHalfWidth + fSelfEnergyTotal;
98 
99  fDeltaPropagator2 = 1.0 / (ReDelta * ReDelta + ImDelta * ImDelta);
100 }
101 
102 /*! <ul>
103  * <li> calculalte p- and s-wave absorption cross section
104  * <li> calculalte total qel (el+cex) p- and s-wave cross section
105  * <li> calculalte fraction of cex in total qel (based on isospin relation)
106  * </ul>
107  */
109 {
110 
111  // common part for all p-wave cross sections
112  const double pXsecCommon = fNormFactor * fCouplingFactor * fDeltaPropagator2 *
113  fMomentumCMS2 / fPionMomentum;
114 
115  // ----- ABSORPTION ----- //
116 
117  // constant factor for p-wave absorption cross section
118  static const double pAborptionFactor = 4.0 / 9.0;
119 
120  // absorption p-wave cross section (see eq. 2.24 and eq. 2.21)
121  const double pXsecAbsorption = pAborptionFactor * pXsecCommon *
122  fSelfEnergyAbsorption;
123 
124  // constant factor for s-wave absorption cross section
125  static const double sAbsorptionFactor = 4.0 * kPi * 197.327 * 10.0 * ImB0;
126 
127  // absorption s-wave cross section (see sec. 3.3)
128  const double sXsecAbsorption = sAbsorptionFactor / fPionMomentum * fNuclearDensity *
129  (1.0 + fPionEnergy / 2.0 / fNucleonMass) / pow (fPionMass / 197.327, 4.0);
130 
131  // total absorption cross section coming from both s- and p-waves
132  fAbsorptionCrossSection = pXsecAbsorption + sXsecAbsorption;
133 
134  // ----- TOTAL ----- //
135 
136  // el+cex p-wave cross section (will be multipled by proper factor later)
137  const double pXsecTotalQel = pXsecCommon * fReducedHalfWidth;
138 
139  // see eq. 3.7
140  const double ksi = (fInvariantMass - fNucleonMass - fPionMass) / fPionMass;
141 
142  // el+cex s-wave cross section
143  const double sXsecTotalQel = quadraticFunction (ksi, fCoefSigma) /
144  fPionMass2 * fNormFactor;
145 
146  // see eq. 3.8
147  const double B = quadraticFunction (ksi, fCoefB);
148  const double D = quadraticFunction (ksi, fCoefD);
149  const double A = 0.5 + 0.5*D;
150  const double C = 1.0 - A;
151 
152  // see 3.4 (chi = 1 for neutron, chi = -1 for proton, chi = 0 for N=Z)
153  const double sTotalQelFactor[fNChannels] = {C - B, A + B, 1.0};
154 
155  // see eq. 2.18 (chi = 1 for neutron, chi = -1 for proton, chi = 0 for N=Z)
156  static const double pTotalQelFactor[fNChannels] = {2.0 / 9.0, 2.0 / 3.0,
157  5.0 / 9.0};
158 
159  // total cross section for each channel = qel_s + qel_p + absorption
160  for (unsigned int i = 0; i < fNChannels; i++)
161  fQelCrossSections[i] = pTotalQelFactor[i] * pXsecTotalQel + sTotalQelFactor[i] * sXsecTotalQel;
162 
163  // ----- CEX ----- //
164 
165  // see 2.18 (chi = 1 for neutron, chi = -1 for proton, chi = 0 for N=Z)
166  static const double pCexFactor[fNChannels] = {4.0 / 27.0, 0.0, 5.0 / 27.0};
167 
168  // see eq. 3.4 (chi = 1 for neutron, chi = -1 for proton, chi = 0 for N=Z)
169  const double sCexFactor[fNChannels] = {2.0 * C, 0.0, 2.0 * C};
170 
171  // total cex cross section coming from both s- and p-waves
172  for (unsigned int i = 0; i < fNChannels; i++)
173  fCexCrossSections[i] = pCexFactor[i] * pXsecTotalQel + sCexFactor[i] * sXsecTotalQel;
174 }
175 
176 //! related to Pauli blocking, see sec. 2.3
178 {
179  // assuming nucleon at rest
180  const double deltaEnergy = fPionEnergy + fNucleonMass;
181  // nucleon energy in CMS
182  const double energyCMS = sqrt(fMomentumCMS * fMomentumCMS + fNucleonMass2);
183  // fermiEnergy calculated from density
184  const double mu0 = (deltaEnergy * energyCMS - fFermiEnergy * fInvariantMass) /
185  fPionMomentum / fMomentumCMS;
186 
187  // see eq. 2.13
188  if (mu0 < -1.0) return 0.0;
189  if (mu0 > 1.0) return 1.0;
190 
191  return (mu0*mu0*mu0 + mu0 + 2) / 4.0;
192 }
193 
194 //! based on eq. 2.21
196 {
197  const double x = fPionKineticEnergy / fPionMass;
198  const double densityFraction = fNuclearDensity / fNormalDensity;
199 
200  const double alpha = quadraticFunction (x, fCoefAlpha);
201  const double beta = quadraticFunction (x, fCoefBeta);
202 
203  double absNN = quadraticFunction (x, fCoefCA2);
204  double absNNN = quadraticFunction (x, fCoefCA3);
205 
206  absNN *= pow (densityFraction, beta);
207 
208  if (absNNN < 0.0) absNNN = 0.0; // it may happen for Tk < 50 MeV
209  else absNNN *= pow (densityFraction, 2.0 * beta);
210 
211  // absNN and absNNN are multiplied by number of nucleons to get proper
212  // cross section normalization
213  //fSelfEnergyAbsorption = 2.0 * absNN + 3.0 * absNNN;
214  fSelfEnergyAbsorption = absNN + absNNN;
215 
216  // this one goes to the delta propagator
217  fSelfEnergyTotal = absNN + absNNN + quadraticFunction (x, fCoefCQ) * pow (densityFraction, alpha);
218 }
219 
220 /*! <ul>
221  * <li> set up nucleus
222  * <li> set up kinematics
223  * <li> set up Delta (width, propagator)
224  * <li> calculate cross sections
225  * </ul>
226  */
227 void INukeOsetFormula :: setupOset (const double &density, const double &pionTk, const int &pionPDG,
228  const double &protonFraction)
229 {
230  setNucleus (density);
231  setKinematics (pionTk, pionPDG == kPdgPi0);
232  setDelta();
233  setCrossSections ();
234  INukeOset::setCrossSections (pionPDG, protonFraction);
235 }
236 
static const double fCoefSigma[fNChannels]
s-wave parametrization eq. 3.7
double deltaReduction() const
calculalte delta width reduction in nuclear medium
void setCrossSections()
calculalte cross sections for each channel
void setKinematics(const double &pionTk, const bool &isPi0)
do kinematics
static const double kNucleonMass
static const double fCoefAlpha[fNChannels]
alpha (eq. 2.21)
static const double fCoefD[fNChannels]
s-wave parametrization eq. 3.8
double quadraticFunction(const double &x, const double *a)
Definition: INukeOset.h:112
static const double fCoefCA2[fNChannels]
two-body absorption (eq. 2.21)
static const double fCoefCQ[fNChannels]
quasi-elastic term (eq. 2.21)
void setDelta()
set up Delta
static const double ImB0
s-wave parametrization eq. 3.12
static const double fNormFactor
MeV^-2 -&gt; mb.
static const double fNormalDensity
normal nuclear density [fm-3]
void setNucleus(const double &density)
set nuclear density and Fermi momentum / energy
static constexpr double A
Definition: Units.h:74
void setSelfEnergy()
calculate delta self energy
static const double fNucleonMass2
average nucleon mass squared
const int kPdgPiP
Definition: PDGCodes.h:158
const int kPdgPi0
Definition: PDGCodes.h:160
static const double fDeltaMass
delta mass in MeV
static const double fCoefCA3[fNChannels]
three-body absorption (eq. 2.21)
static const double fNucleonMass
average nucleon mass [MeV]
static const double fCoefB[fNChannels]
s-wave parametrization eq. 3.8
virtual void setCrossSections()=0
calculalte cross sections for each channel
void setupOset(const double &density, const double &pionTk, const int &pionPDG, const double &protonFraction)
use to set up Oset class (assign pion Tk, nuclear density etc)
static const double fCouplingConstant
f*^2
static const double fCoefBeta[fNChannels]
beta (eq. 2.21)