14 #include "Framework/Conventions/GBuild.h"
25 using namespace genie;
26 using namespace genie::constants;
27 using namespace genie::utils;
57 const Kinematics & kinematics = interaction -> Kine();
58 const InitialState & init_state = interaction -> InitState();
64 double Q2 = kinematics.
Q2();
65 double y = kinematics.
y();
66 double x = kinematics.
x();
72 LOG(
"BergerSehgalCohPi",
pDEBUG) <<
"Pion COM momentum negative for Q2 = " << Q2 <<
73 " x = " << x <<
" y = " << y;
78 LOG(
"BergerSehgalCohPi",
pDEBUG) <<
"Exact kin. form = " << front <<
79 " E = " << E <<
" Q2 = " << Q2 <<
" y = " << y <<
" x = " << x;
83 double A = (double) init_state.
Tgt().
A();
84 double A2 = TMath::Power(A, 2.);
85 double A_3 = TMath::Power(A, 1./3.);
86 double M = init_state.
Tgt().
Mass();
89 double ma2 = TMath::Power(
fMa, 2);
90 double Ga = ma2 / (ma2 +
Q2);
91 double Ga2 = TMath::Power(Ga, 2.);
96 double Epi2 = TMath::Power(Epi, 2.);
98 double R2 = TMath::Power(R, 2.);
99 double b = 0.33333 * R2;
100 double MxEpi = M * x / Epi;
101 double mEpi2 = (M_pi * M_pi) / Epi2;
102 double tA = 1. + MxEpi - 0.5 * mEpi2;
103 double tB = TMath::Sqrt(1.0 + 2 * MxEpi) * TMath::Sqrt(1.0 - mEpi2);
104 double tmin = 2 * Epi2 * (tA - tB);
105 double tmax = 2 * Epi2 * (tA + tB);
115 double siginel_pin = sigtot_pin - sigel_pin;
120 double fabs_input = (9.0 * A_3) / (16.0 *
kPi * Ro2);
121 double fabs = TMath::Exp( -1.0 * fabs_input * siginel_pin);
128 double RS_factor = (A2 * fabs) / (16.0 *
kPi) * (sigtot_pin * sigtot_pin);
131 double tpi = (E * y) - M_pi - ((Q2 + M_pi * M_pi) / (2 * M));
134 double tpihigh = 0.0;
135 double sighigh = 0.0;
136 double dsigdzfit = 0.0;
137 double dsigdtfit = 0.0;
141 double logt_step = TMath::Abs(log(tmax) - log(tmin)) / tstep;
142 double logt = log(tmin) - logt_step/2.0;
143 double t_itt = TMath::Exp(logt);
144 double t_width = 0.0;
146 for (
double t_step = 0; t_step<tstep; t_step++) {
148 logt = logt + logt_step;
149 t_itt = TMath::Exp(logt);
150 t_width = t_itt*logt_step;
155 LOG(
"BergerSehgalCohPi",
pERROR) <<
"Call to PionNucleusXSec code failed - return xsec of 0.0";
158 dsigdzfit = siglow + (sighigh - siglow) * (tpi - tpilow) / (tpihigh - tpilow);
159 dsigdtfit = dsigdzfit / (2.0 * ppistar * ppistar);
161 dsig += 1.0 * front * Ga2 * t_width * dsigdtfit *
units::mb;
164 dsig += front * Ga2 * t_width * RS_factor * exp(-1.0*b*t_itt);
179 double ml2 = TMath::Power(ml,2);
180 double Q2min = ml2 * y/(1-y);
182 double C1 = TMath::Power(Ga - 0.5 * Q2min / (Q2 +
kPionMass2), 2);
183 double C2 = 0.25 * y * Q2min * (Q2 - Q2min) /
192 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
194 <<
"\n momentum transfer .............. Q2 = " << Q2
195 <<
"\n mass number .................... A = " << A
196 <<
"\n pion energy .................... Epi = " << Epi
197 <<
"\n propagator term ................ Ga2 = " << Ga2
198 <<
"\n total pi+N cross section ....... sigT = " << sigtot_pin
199 <<
"\n inelastic pi+N cross section ... sigI = " << siginel_pin
200 <<
"\n nuclear size scale ............. Ro = " <<
fRo
201 <<
"\n pion absorption factor ......... Fabs = " << fabs
202 <<
"\n t integration range ............ [" << tmin <<
"," << tmax <<
"]"
204 <<
"d2xsec/dQ2dy[COHPi] (Q2= " << Q2 <<
", y="
205 << y <<
", E=" << E <<
") = "<< xsec;
222 const Kinematics & kinematics = interaction -> Kine();
223 const InitialState & init_state = interaction -> InitState();
228 double Q2 = kinematics.
Q2();
229 double y = kinematics.
y();
230 double fp2 = (0.93 * M_pi)*(0.93 * M_pi);
232 double term = ((
kGF2 * fp2) / (4.0 *
kPi2)) *
233 ((E * (1.0 - y)) / sqrt(y*E * y*E + Q2)) *
234 (1.0 - Q2 / (4.0 * E*E * (1.0 - y)));
242 const Kinematics & kinematics = interaction -> Kine();
243 const InitialState & init_state = interaction -> InitState();
248 double Q2 = kinematics.
Q2();
249 double y = kinematics.
y();
250 double MT = init_state.
Tgt().
Mass();
252 double W2 = MT*MT - Q2 + 2.0 * y * E * MT;
253 double arg = (2.0*MT*(y*E - M_pi) - Q2 - M_pi*M_pi)*(2.0*MT*(y*E + M_pi) - Q2 - M_pi*M_pi);
254 if (arg < 0)
return arg;
255 double ppistar = TMath::Sqrt(arg) / 2.0 / TMath::Sqrt(W2);
272 const Target & target = init_state.
Tgt();
277 if (!proc_info.
IsWeak())
return false;
279 if (!(target.
A()>1))
return false;
304 fCos8c2 = TMath::Power(TMath::Cos(thc), 2);
Cross Section Calculation Interface.
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
double J(double q0, double q3, double Enu, double ml)
Cross Section Integrator Interface.
bool fRSPionXSec
Use Rein-Sehgal "style" pion-nucleon xsecs.
double Q2(const Interaction *const i)
static const double kPi0Mass
Generated/set kinematical variables for an event.
double PionNucleonXSec(double Epion, bool get_total, bool isChargedPion=true)
double fRo
nuclear size scale parameter
virtual ~BergerSehgalCOHPiPXSec2015()
double x(bool selected=false) const
bool IsCoherentProduction(void) const
enum genie::EKinePhaseSpace KinePhaseSpace_t
double PionCOMAbsMomentum(const Interaction *i) const
double y(bool selected=false) const
static constexpr double b
BergerSehgalCOHPiPXSec2015()
Summary information for an interaction.
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
static constexpr double A
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
static constexpr double mb
bool IsAntiNeutrino(int pdgc)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
virtual void Configure(const Registry &config)
double ExactKinematicTerm(const Interaction *i) const
double fCos8c2
cos^2(Cabibbo angle)
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
void Configure(const Registry &config)
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const XSecIntegratorI * fXSecIntegrator
static const double kPionMass
bool HitNucIsSet(void) const
A registry. Provides the container for algorithm configuration parameters.
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
static constexpr double fermi
const InitialState & InitState(void) const
const ProcessInfo & ProcInfo(void) const
double Q2(bool selected=false) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
const UInt_t kISkipProcessChk
if set, skip process validity checks
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double Integral(const Interaction *i) const
Initial State information.
int PionNucleusXSec(double tpi, double ppistar, double t_new, double A, double &tpilow, double &siglow, double &tpihigh, double &sighigh)
static const double kPionMass2
const Algorithm * SubAlg(const RgKey ®istry_key) const