29 using namespace genie;
30 using namespace genie::constants;
57 const InitialState & init_state = interaction -> InitState();
58 const Kinematics & kinematics = interaction -> Kine();
62 double E2 = TMath::Power(E,2);
63 double Q2 = kinematics.
Q2();
64 double W = kinematics.
W();
66 double MN2 = TMath::Power(MN,2);
70 LOG(
"PaschLal",
pDEBUG) <<
"Input kinematics: W = " << W <<
", Q2 = " <<
Q2;
75 double MR2 = TMath::Power(MR,2);
76 double MR3 = TMath::Power(MR,3);
80 const double kPlRes_f3_P1232_V = 1.95/MN;
81 const double kPlRes_f4_P1232_V = -1.95/MN;
82 const double kPlRes_f5_P1232_V = 0;
83 const double kPlRes_f5_P1232_A = 1.2;
84 const double kPlRes_f4_P1232_A = -0.3/MN2;
85 const double kPlRes_f3_P1232_A = 0;
86 const double kPlRes_f6_P1232_A = kPlRes_f5_P1232_A;
88 double MA2 = TMath::Power(
fMa, 2 );
89 double MV2 = TMath::Power(
fMv, 2 );
90 double ftmp1a = TMath::Power( 1 + Q2/MA2, 2 );
91 double ftmp1v = TMath::Power( 1 + Q2/MA2, 2 );
92 double ftmp2a = 1 + Q2/3/MA2;
93 double ftmp2v = 1 + Q2/4/MV2;
94 double f3A = kPlRes_f3_P1232_A/ftmp1a/ftmp2a;
95 double f4A = kPlRes_f4_P1232_A/ftmp1a/ftmp2a;
96 double f5A = kPlRes_f5_P1232_A/ftmp1a/ftmp2a;
97 double f6A = kPlRes_f6_P1232_A/ftmp1a/ftmp2a/(Q2+Mpi2);
98 double f3V = kPlRes_f3_P1232_V/ftmp1v/ftmp2v;
99 double f4V = kPlRes_f4_P1232_V/ftmp1v/ftmp2v/
W;
100 double f5V = kPlRes_f5_P1232_V/ftmp1v/ftmp2v;
101 double f3V4A = f3V*f4A;
102 double f3V5A = f3V*f5A;
103 double f4V4A = f4V*f4A;
104 double f4V5A = f4V*f5A;
105 double f3A2 = TMath::Power( f3A, 2 );
106 double f4A2 = TMath::Power( f4A, 2 );
107 double f5A2 = TMath::Power( f5A, 2 );
108 double f6A2 = TMath::Power( f6A, 2 );
109 double f3V2 = TMath::Power( f3V, 2 );
110 double f4V2 = TMath::Power( f4V, 2 );
111 double f5V2 = TMath::Power( f5V, 2 );
117 double Gamma_R=Gamma_R0*pow((this->
PPiStar(W,MN)/this->
PPiStar(MR,MN)),3);
120 if (
GetConfig().Exists(
"running-gamma") ) {
124 if ( gamma_model.find(
"Hagiwara") != string::npos )
127 Gamma_R = Gamma_R0*MR/W*pow((this->
PPiStar(W,MN)/this->
PPiStar(MR,MN)),1);
129 if ( gamma_model.find(
"Galster") != string::npos )
132 double gtmp1 = TMath::Power( this->
PPiStar(W,MN) / this->
PPiStar(MR,MN), 3);
133 double gtmp2 = TMath::Power( this->
PPiStar(W,MN) / this->
PPiStar(MR,MN), 2);
134 Gamma_R = Gamma_R0*gtmp1/(1+gtmp2);
137 double Breit_Wigner = TMath::Power(W*W-MR2,2) + MR2 * TMath::Power(Gamma_R,2);
144 double nu = this->
Nu(Q2,W,MN);
146 double qk = -(Q2+Mmu2)/2.;
148 double pq2 = TMath::Power(pq,2);
149 double pq3 = TMath::Power(pq,3);
150 double Q4 = TMath::Power(Q2,2);
153 double W1=0, W2=0, W3=0, W4=0, W5=0;
155 W1 = 3.*(2*f5A2*MN2*MR2+2*f5A2*MN*MR3+2*f3A*f5A*MN2*MR*pq+2*f5A2*MR2*pq
156 +4*f3A*f5A*MN*MR2*pq+4*f4A*f5A*MN2*MR2*pq+2*f3A*f5A*MR3*pq
157 +4*f4A*f5A*MN*MR3*pq+2*f3A2*MN2*pq2+2*f3V2*MN2*pq2+2*f3A*f5A*MR*pq2
158 +2*f3A*f4A*MN2*MR*pq2+2*f3V*f4V*MN2*MR*pq2+2*f3V*f5V*MN2*MR*pq2
159 +2*f3A2*MR2*pq2+2*f3V2*MR2*pq2+4*f4A*f5A*MR2*pq2+4*f3A*f4A*MN*MR2*pq2
160 -4*f3V*f4V*MN*MR2*pq2-4*f3V*f5V*MN*MR2*pq2+2*f4A2*MN2*MR2*pq2+2*f4V2*MN2*MR2*pq2
161 +4*f4V*f5V*MN2*MR2*pq2+2*f5V2*MN2*MR2*pq2+2*f3A*f4A*MR3*pq2+2*f3V*f4V*MR3*pq2
162 +2*f3V*f5V*MR3*pq2+2*f4A2*MN*MR3*pq2-2*f4V2*MN*MR3*pq2-4*f4V*f5V*MN*MR3*pq2
163 -2*f5V2*MN*MR3*pq2+2*f3A2*pq3+2*f3V2*pq3+2*f3A*f4A*MR*pq3+2*f3V*f4V*MR*pq3
164 +2*f3V*f5V*MR*pq3+2*f4A2*MR2*pq3+2*f4V2*MR2*pq3+4*f4V*f5V*MR2*pq3+2*f5V2*MR2*pq3
165 -2*f3A*f5A*MN2*MR*Q2-4*f3A*f5A*MN*MR2*Q2+2*f3A2*MN2*MR2*Q2
166 +2*f3V2*MN2*MR2*Q2-4*f4A*f5A*MN2*MR2*Q2-2*f3A2*MN*MR3*Q2+2*f3V2*MN*MR3*Q2
167 -4*f4A*f5A*MN*MR3*Q2-4*f3A2*MN2*pq*Q2-4*f3V2*MN2*pq*Q2-2*f3A*f5A*MR*pq*Q2
168 -4*f3A*f4A*MN2*MR*pq*Q2-4*f3V*f4V*MN2*MR*pq*Q2-2*f3V*f5V*MN2*MR*pq*Q2
169 -4*f4A*f5A*MR2*pq*Q2-8*f3A*f4A*MN*MR2*pq*Q2+8*f3V*f4V*MN*MR2*pq*Q2
170 +4*f3V*f5V*MN*MR2*pq*Q2-4*f4A2*MN2*MR2*pq*Q2-4*f4V2*MN2*MR2*pq*Q2
171 -4*f4V*f5V*MN2*MR2*pq*Q2-2*f3A*f4A*MR3*pq*Q2-2*f3V*f4V*MR3*pq*Q2
172 -4*f4A2*MN*MR3*pq*Q2+4*f4V2*MN*MR3*pq*Q2+4*f4V*f5V*MN*MR3*pq*Q2
173 -4*f3A2*pq2*Q2-4*f3V2*pq2*Q2-4*f3A*f4A*MR*pq2*Q2-4*f3V*f4V*MR*pq2*Q2
174 -2*f3V*f5V*MR*pq2*Q2-4*f4A2*MR2*pq2*Q2-4*f4V2*MR2*pq2*Q2-4*f4V*f5V*MR2*pq2*Q2
175 +2*f3A2*MN2*Q4+2*f3V2*MN2*Q4+2*f3A*f4A*MN2*MR*Q4+2*f3V*f4V*MN2*MR*Q4
176 +4*f3A*f4A*MN*MR2*Q4-4*f3V*f4V*MN*MR2*Q4+2*f4A2*MN2*MR2*Q4+2*f4V2*MN2*MR2*Q4
177 +2*f4A2*MN*MR3*Q4-2*f4V2*MN*MR3*Q4+2*f3A2*pq*Q4+2*f3V2*pq*Q4+2*f3A*f4A*MR*pq*Q4
178 +2*f3V*f4V*MR*pq*Q4+2*f4A2*MR2*pq*Q4+2*f4V2*MR2*pq*Q4)/(3.*MR2);
181 +f5A2*MN*MR+f5A2*pq+f3A2*MN2*Q2+f3V2*MN2*Q2+f3A*f5A*MR*Q2+f3A*f4A*MN2*MR*Q2
182 +f3V*f4V*MN2*MR*Q2+f3V*f5V*MN2*MR*Q2+f3A2*MR2*Q2+f3V2*MR2*Q2
183 +2*f3A*f4A*MN*MR2*Q2-2*f3V*f4V*MN*MR2*Q2-2*f3V*f5V*MN*MR2*Q2
184 +f4A2*MN2*MR2*Q2+f4V2*MN2*MR2*Q2+2*f4V*f5V*MN2*MR2*Q2+f5V2*MN2*MR2*Q2+f3A*f4A*MR3*Q2
185 +f3V*f4V*MR3*Q2+f3V*f5V*MR3*Q2+f4A2*MN*MR3*Q2-f4V2*MN*MR3*Q2
186 -2*f4V*f5V*MN*MR3*Q2-f5V2*MN*MR3*Q2+f3A2*pq*Q2+f3V2*pq*Q2+f3A*f4A*MR*pq*Q2
187 +f3V*f4V*MR*pq*Q2+f3V*f5V*MR*pq*Q2+f4A2*MR2*pq*Q2+f4V2*MR2*pq*Q2
188 +2*f4V*f5V*MR2*pq*Q2+f5V2*MR2*pq*Q2+f5V2*MN2*Q4+f3V*f5V*MR*Q4
189 -f5V2*MN*MR*Q4+f5V2*pq*Q4))/(3.*MR2);
191 W3 = 3.*((f3V4A*(Q2-pq)-f3V5A)*(2*MR2+2*MN*MR+Q2-pq)*4./3./MR
192 -(Q2-pq)*(f4V4A*(Q2-pq)-f4V5A)*4./3.);
195 W4 = 3.*(2*(f5A2*MN2+f5A2*MN*MR+f3A*f5A*MN2*MR
196 +2*f3A*f5A*MN*MR2-f3A2*MN2*MR2-f3V2*MN2*MR2+2*f4A*f5A*MN2*MR2-2*f5A*f6A*MN2*MR2
197 +f3A2*MN*MR3-f3V2*MN*MR3+2*f4A*f5A*MN*MR3-2*f5A*f6A*MN*MR3+f5A2*pq
198 +2*f3A2*MN2*pq+2*f3V2*MN2*pq+2*f5A*f6A*MN2*pq+2*f3A*f5A*MR*pq
199 +2*f5A*f6A*MN*MR*pq+2*f3A*f4A*MN2*MR*pq+2*f3V*f4V*MN2*MR*pq
200 +f3V*f5V*MN2*MR*pq-f3A*f6A*MN2*MR*pq+2*f4A*f5A*MR2*pq
201 -2*f5A*f6A*MR2*pq+4*f3A*f4A*MN*MR2*pq-4*f3V*f4V*MN*MR2*pq
202 -2*f3V*f5V*MN*MR2*pq-2*f3A*f6A*MN*MR2*pq+2*f4A2*MN2*MR2*pq
203 +2*f4V2*MN2*MR2*pq+2*f4V*f5V*MN2*MR2*pq-2*f4A*f6A*MN2*MR2*pq+f3A*f4A*MR3*pq
204 +f3V*f4V*MR3*pq-f3A*f6A*MR3*pq+2*f4A2*MN*MR3*pq-2*f4V2*MN*MR3*pq
205 -2*f4V*f5V*MN*MR3*pq-2*f4A*f6A*MN*MR3*pq+2*f3A2*pq2+2*f3V2*pq2
206 +2*f5A*f6A*pq2+f5V2*MN2*pq2+f6A2*MN2*pq2+2*f3A*f4A*MR*pq2+2*f3V*f4V*MR*pq2
207 +2*f3V*f5V*MR*pq2-f5V2*MN*MR*pq2+f6A2*MN*MR*pq2+2*f4A2*MR2*pq2+2*f4V2*MR2*pq2
208 +2*f4V*f5V*MR2*pq2-2*f4A*f6A*MR2*pq2+f5V2*pq3+f6A2*pq3-f3A2*MN2*Q2-f3V2*MN2*Q2
209 -2*f5A*f6A*MN2*Q2-2*f5A*f6A*MN*MR*Q2-f3A*f4A*MN2*MR*Q2
210 -f3V*f4V*MN2*MR*Q2-2*f3A*f4A*MN*MR2*Q2+2*f3V*f4V*MN*MR2*Q2
211 -f4A2*MN2*MR2*Q2-f4V2*MN2*MR2*Q2+f6A2*MN2*MR2*Q2-f4A2*MN*MR3*Q2+f4V2*MN*MR3*Q2
212 +f6A2*MN*MR3*Q2-f3A2*pq*Q2-f3V2*pq*Q2-2*f5A*f6A*pq*Q2-2*f6A2*MN2*pq*Q2
213 -f3A*f4A*MR*pq*Q2-f3V*f4V*MR*pq*Q2-f3A*f6A*MR*pq*Q2
214 -2*f6A2*MN*MR*pq*Q2-f4A2*MR2*pq*Q2-f4V2*MR2*pq*Q2+f6A2*MR2*pq*Q2
215 -2*f6A2*pq2*Q2+f6A2*MN2*Q4+f6A2*MN*MR*Q4+f6A2*pq*Q4))/(3.*MR2);
218 +2*f5A2*MN*MR+f3A*f5A*MN2*MR+2*f3A*f5A*MN*MR2+2*f4A*f5A*MN2*MR2
219 +f3A*f5A*MR3+2*f4A*f5A*MN*MR3+2*f5A2*pq+2*f3A2*MN2*pq+2*f3V2*MN2*pq
220 +2*f5A*f6A*MN2*pq+2*f3A*f5A*MR*pq+2*f5A*f6A*MN*MR*pq
221 +2*f3A*f4A*MN2*MR*pq+2*f3V*f4V*MN2*MR*pq+2*f3V*f5V*MN2*MR*pq
222 +2*f3A2*MR2*pq+2*f3V2*MR2*pq+2*f4A*f5A*MR2*pq+4*f3A*f4A*MN*MR2*pq
223 -4*f3V*f4V*MN*MR2*pq-4*f3V*f5V*MN*MR2*pq+2*f4A2*MN2*MR2*pq
224 +2*f4V2*MN2*MR2*pq+4*f4V*f5V*MN2*MR2*pq+2*f5V2*MN2*MR2*pq+2*f3A*f4A*MR3*pq
225 +2*f3V*f4V*MR3*pq+2*f3V*f5V*MR3*pq+2*f4A2*MN*MR3*pq-2*f4V2*MN*MR3*pq
226 -4*f4V*f5V*MN*MR3*pq-2*f5V2*MN*MR3*pq+2*f3A2*pq2+2*f3V2*pq2+2*f5A*f6A*pq2
227 +2*f3A*f4A*MR*pq2+2*f3V*f4V*MR*pq2+2*f3V*f5V*MR*pq2+2*f4A2*MR2*pq2
228 +2*f4V2*MR2*pq2+4*f4V*f5V*MR2*pq2+2*f5V2*MR2*pq2-2*f5A*f6A*MN2*Q2+f3A*f5A*MR*Q2
229 -2*f5A*f6A*MN*MR*Q2-f3A*f6A*MN2*MR*Q2-2*f3A*f6A*MN*MR2*Q2
230 -2*f4A*f6A*MN2*MR2*Q2-f3A*f6A*MR3*Q2-2*f4A*f6A*MN*MR3*Q2
231 -2*f5A*f6A*pq*Q2+2*f5V2*MN2*pq*Q2+2*f3V*f5V*MR*pq*Q2
232 -2*f5V2*MN*MR*pq*Q2-2*f4A*f6A*MR2*pq*Q2+2*f5V2*pq2*Q2-f3A*f6A*MR*Q4)/(3.*MR2);
234 double s1 = W1 * (Q2+Mmu2)
235 + W2 * (2*pk*pk-2*pq*pk+MN*qk)
237 + W4 * Mmu2*(Q2+Mmu2)/2.
254 int NNucl = (isp) ? target.
Z() : target.
N();
289 double fermi_constant ;
290 GetParam(
"FermiConstant", fermi_constant ) ;
298 fCos28c = TMath::Power( TMath::Cos(thc), 2 );
317 double p_pi_star = this->
PPiStar(W,MN);
318 double nu_star = this->
NuStar(Q2,W,MN);
320 double p_pi_star_2 = TMath::Power(p_pi_star, 2);
321 double p_pi_star_4 = TMath::Power(p_pi_star_2, 2);
322 double nu_star_2 = TMath::Power(nu_star, 2);
324 double q3 = TMath::Sqrt(Q2+nu_star_2);
325 double q6 = TMath::Power(q3,2);
326 double q12 = TMath::Power(q6,2);
328 double qF2 = TMath::Power(qF,2);
329 double qF3 = TMath::Power(qF,3);
331 if ( q3+p_pi_star < 2*qF )
333 Paulii = ( (3*q6 + p_pi_star_2)/2/qF
334 -(5*q12+p_pi_star_4 + 10*q6*p_pi_star_2)/40/qF3
337 if ( (q3+p_pi_star > 2*qF) && (q3-p_pi_star < 2*qF) )
339 double tmp1 = TMath::Power( q3+p_pi_star, 2 );
340 double tmp2 = TMath::Power( q3-p_pi_star, 3 );
341 double tmp3 = TMath::Power( q3-p_pi_star, 5 );
343 Paulii = (tmp1-4.0*qF2/5.0 - tmp2/2/qF + tmp3/40/qF3)/4/p_pi_star/q3;
345 if ( q3-p_pi_star > 2*qF )
355 return (TMath::Power(W,2) - TMath::Power(MN,2) + Q2)/2/MN;
360 double W2 = TMath::Power(W,2);
364 return TMath::Sqrt( (W2-a)*(W2-b) )/2/
W;
369 return (TMath::Power(W,2) - TMath::Power(MN,2) - Q2)/2/
W;
bool fTurnOnPauliCorrection
P33PaschosLalakulichPXSec()
Cross Section Calculation Interface.
double W(bool selected=false) const
double J(double q0, double q3, double Enu, double ml)
double PPiStar(double W, double MN) const
...
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
int HitNucPdg(void) const
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double HitNucMass(void) const
Generated/set kinematical variables for an event.
void Configure(const Registry &config)
double Mass(Resonance_t res)
resonance mass (GeV)
double NuStar(double Q2, double W, double MN) const
...
double Width(Resonance_t res)
resonance width (GeV)
enum genie::EKinePhaseSpace KinePhaseSpace_t
virtual const Registry & GetConfig(void) const
double Nu(double Q2, double W, double MN) const
kinematic variables
static constexpr double b
double W(const Interaction *const i)
double Pauli(double Q2, double W, double MN) const
Pauli suppression for D2.
virtual ~P33PaschosLalakulichPXSec()
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...
const XSecIntegratorI * fXSecIntegrator
static const double kMuonMass2
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
virtual void Configure(const Registry &config)
double Integral(const Interaction *i) const
RgStr GetString(RgKey key) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
static const double kPionMass
A registry. Provides the container for algorithm configuration parameters.
const UInt_t kIAssumeFreeNucleon
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
double Q2(bool selected=false) const
bool GetParamDef(const RgKey &name, T &p, const T &def) 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
const UInt_t kISkipProcessChk
if set, skip process validity checks
Initial State information.
static const double kPionMass2
const Algorithm * SubAlg(const RgKey ®istry_key) const