GENIEGenerator
Loading...
Searching...
No Matches
MKSPPPXSec2020.cxx
Go to the documentation of this file.
1//____________________________________________________________________________
2/*
3 Copyright (c) 2003-2025, The GENIE Collaboration
4 For the full text of the license visit http://copyright.genie-mc.org
5 or see $GENIE/LICENSE
6
7 Authors: Igor Kakorin <kakorin@jinr.ru>, Joint Institute for Nuclear Research
8 Vadim Naumov <vnaumov@theor.jinr.ru >, Joint Institute for Nuclear Research \n
9 adapted from code provided by
10 Minoo Kabirnezhad <minoo.kabirnezhad@physics.ox.ac.uk>
11 University of Oxford, Department of Physics
12 based on code of Costas Andreopoulos <c.andreopoulos \at cern.ch>
13 University of Liverpool
14
15
16 For the class documentation see the corresponding header file.
17*/
18//____________________________________________________________________________
19
20#include <TMath.h>
21#include <TSystem.h>
22
23#include "Framework/Algorithm/AlgFactory.h"
24#include "Framework/Algorithm/AlgConfigPool.h"
25#include "Framework/ParticleData/BaryonResUtils.h"
26#include "Framework/Conventions/GBuild.h"
27#include "Framework/Conventions/Constants.h"
28#include "Framework/Conventions/RefFrame.h"
29#include "Framework/Conventions/KineVar.h"
30#include "Framework/Conventions/Units.h"
31#include <Framework/Conventions/KinePhaseSpace.h>
32#include "Framework/Messenger/Messenger.h"
33#include "Framework/ParticleData/PDGCodes.h"
34#include "Framework/ParticleData/PDGUtils.h"
35#include "Framework/ParticleData/PDGLibrary.h"
36#include "Framework/Interaction/SppChannel.h"
37#include "Framework/Utils/KineUtils.h"
38#include "Framework/Utils/Range1.h"
39#include "Physics/XSectionIntegration/XSecIntegratorI.h"
40#include "Physics/Resonance/XSection/MKSPPPXSec2020.h"
41#include "Physics/NuclearState/FermiMomentumTablePool.h"
42#include "Physics/NuclearState/FermiMomentumTable.h"
43#include "Physics/NuclearState/NuclearUtils.h"
44
45#include <algorithm>
46
47
48using namespace genie;
49using namespace genie::constants;
50
51//____________________________________________________________________________
52MKSPPPXSec2020::MKSPPPXSec2020() :
53XSecAlgorithmI("genie::MKSPPPXSec2020")
54{
55
56}
57//____________________________________________________________________________
58MKSPPPXSec2020::MKSPPPXSec2020(string config) :
59XSecAlgorithmI("genie::MKSPPPXSec2020", config)
60{
61
62}
63//____________________________________________________________________________
68//____________________________________________________________________________
69double MKSPPPXSec2020::XSec(const Interaction * interaction, KinePhaseSpace_t kps) const
70{
71 if(! this -> ValidProcess (interaction) ) return 0;
72 if(! this -> ValidKinematics (interaction) ) return 0;
73
74 const InitialState & init_state = interaction -> InitState();
75 const ProcessInfo & proc_info = interaction -> ProcInfo();
76 const Target & target = init_state.Tgt();
77
78 double xsec = 0;
79 //-- Get 1pi exclusive channel
80 SppChannel_t spp_channel = SppChannel::FromInteraction(interaction);
81
82 int nucpdgc = target.HitNucPdg();
83 int probepdgc = init_state.ProbePdg();
84 bool is_nu = pdg::IsNeutrino (probepdgc);
85 bool is_nubar = pdg::IsAntiNeutrino (probepdgc);
86 bool is_CC = proc_info.IsWeakCC();
87 bool is_NC = proc_info.IsWeakNC();
88 bool is_p = pdg::IsProton (nucpdgc);
89
90
91 // Get kinematical parameters
92 const Kinematics & kinematics = interaction -> Kine();
93 double E = init_state.ProbeE(kRfHitNucRest);
94 double ml = interaction->FSPrimLepton()->Mass();
95 double ml2 = ml*ml;
96 double Q2 = kinematics.Q2();
97 double W = kinematics.W();
98 double W2 = W*W;
99 double Wt2 = W*2;
100
101
102 // dimension of kine phase space
103 std::string s = KinePhaseSpace::AsString(kps);
104 int kpsdim = s!="<|E>"?1 + std::count(s.begin(), s.begin()+s.find('}'), ','):0;
105 if (kpsdim < 3 || kpsdim > 4) return 0;
106 // Pion angles should be given in Adler frame
107 double CosTheta = kinematics.GetKV(kKVctp);
108 double Phi = 0;
109 if (kpsdim == 4)
110 Phi = kinematics.GetKV(kKVphip);
111
112 double SinTheta = TMath::Sqrt(1 - CosTheta*CosTheta);
113 double SinHalfTheta = TMath::Sqrt((1 - CosTheta)/2);
114 double CosHalfTheta = TMath::Sqrt((1 + CosTheta)/2);
115
116 PDGLibrary * pdglib = PDGLibrary::Instance();
117
118 // imply isospin symmetry
119 double m_pi = (pdglib->Find(kPdgPiP)->Mass() + pdglib->Find(kPdgPi0)->Mass() + pdglib->Find(kPdgPiM)->Mass())/3;
120 double m_pi2 = m_pi*m_pi;
121
122 double M = (pdglib->Find(kPdgProton)->Mass() + pdglib->Find(kPdgNeutron)->Mass())/2;
123 double M2 = M*M;
124 double Mt2 = M*2;
125
126 double q_0 = (W2 - M2 + m_pi2)/Wt2;
127 double q_02 = q_0*q_0;
128 double k_0 = (W2 - M2 - Q2)/Wt2;
129 double abs_mom_q = TMath::Sqrt(TMath::Max(0., q_02 - m_pi2));
130 double abs_mom_k = TMath::Sqrt(k_0*k_0 + Q2);
131 //double E_2L = (M2 - W2 - Q2 + Mt2*E)/Mt2;
132 double abs_mom_k_L = W*abs_mom_k/M;
133 double abs_mom_k_L2 = abs_mom_k_L*abs_mom_k_L;
134 double k_2 = (M2 + Mt2*E - W2 - ml2)/Wt2;
135 double k_1 = k_2 + k_0;
136 double p_10 = (W2 + M2 + Q2)/Wt2;
137 double qk = q_0*k_0 - abs_mom_k*abs_mom_q*CosTheta;
138 //double k_2L = TMath::Sqrt(E_2L*E_2L - ml2); //magnitude of lepton momentum in lab frame
139
140 // There is no check of positivity under root expression in the original code
141 double k_2_iso = TMath::Sqrt(TMath::Max(0., k_2*k_2 - ml2)); //magnitude of lepton momentum in isobaric frame
142 double cos_theta = k_2_iso>0?(2*k_1*k_2 - Q2 - ml2)/2/k_1/k_2_iso:0;
143 // There are no checks presented below in the original code
144 if (cos_theta > 1)
145 cos_theta = 1;
146 if (cos_theta < -1)
147 cos_theta = -1;
148
149 // Eq. 7 of ref. 1
150 double A_plus = TMath::Sqrt( k_1*(k_2 - k_2_iso) );
151 double A_minus = TMath::Sqrt( k_1*(k_2 + k_2_iso) );
152 //Eq. 6 of ref. 1
153 double eps_1_plus = 2*A_plus *(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
154 double eps_1_minus = -2*A_minus*(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
155 double eps_2_plus = 2*A_plus *TMath::Sqrt(1 + cos_theta);
156 double eps_2_minus = 2*A_minus*TMath::Sqrt(1 - cos_theta);
157 //Eq. 9 of ref. 1
158 double eps_zero_L = -2*A_minus*TMath::Sqrt(1 + cos_theta); // L->lambda = -1
159 double eps_zero_R = 2*A_plus *TMath::Sqrt(1 - cos_theta); // R->lambda = +1
160 double eps_z_L = -2*A_minus*(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
161 double eps_z_R = 2*A_plus *(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
162
163 if (is_nubar)
164 {
165 if (fUseAuthorCode)
166 {
167 //This ``recipe'' of transition from neutrino to antineutrino case is promoted by Minoo Kabirnezhad.
168 //However, it is not correct in our opinion. All details can be found in Ref. [12], see
169 //section "Problem with transition from neutrino to antineutrino case".
170 Phi = -Phi;
171 eps_1_plus = -2*A_minus *(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
172 eps_1_minus = -2*A_plus*(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
173 eps_2_plus = -2*A_minus *TMath::Sqrt(1 - cos_theta);
174 eps_2_minus = 2*A_plus*TMath::Sqrt(1 + cos_theta);
175 eps_zero_L = 2*A_plus*TMath::Sqrt(1 - cos_theta); // L->lambda = -1
176 eps_zero_R = 2*A_minus *TMath::Sqrt(1 + cos_theta); // R->lambda = +1
177 eps_z_L = 2*A_plus*(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
178 eps_z_R = 2*A_minus *(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
179 }
180 else
181 {
182 eps_1_plus = 2*A_minus *(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
183 eps_1_minus = 2*A_plus*(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
184 eps_2_plus = 2*A_minus *TMath::Sqrt(1 - cos_theta);
185 eps_2_minus = -2*A_plus*TMath::Sqrt(1 + cos_theta);
186 eps_zero_L = 2*A_plus*TMath::Sqrt(1 - cos_theta); // L->lambda = -1
187 eps_zero_R = 2*A_minus *TMath::Sqrt(1 + cos_theta); // R->lambda = +1
188 eps_z_L = 2*A_plus*(k_1 + k_2_iso)/abs_mom_k*TMath::Sqrt(1 - cos_theta);
189 eps_z_R = 2*A_minus *(k_1 - k_2_iso)/abs_mom_k*TMath::Sqrt(1 + cos_theta);
190 }
191 }
192 //Eq. 10 of ref. 1
193 double C_L_plus = k1_Sqrt2*(eps_1_plus - eps_2_plus);
194 double C_L_minus = k1_Sqrt2*(eps_1_minus - eps_2_minus);
195 double C_R_plus = -k1_Sqrt2*(eps_1_plus + eps_2_plus);
196 double C_R_minus = -k1_Sqrt2*(eps_1_minus + eps_2_minus);
197 double C_S_plus = TMath::Sqrt(TMath::Abs(eps_zero_R*eps_zero_R - eps_z_R*eps_z_R));
198 double C_S_minus = TMath::Sqrt(TMath::Abs(eps_zero_L*eps_zero_L - eps_z_L*eps_z_L));
199
200 // if C_S_plus or C_S_minus are equal to zero, then the singularity appers
201 // in the expressions for Hbkg and Hres at BosonPolarization=PLUS0 or BosonPolarization=MINUS0
202 // which further eliminated by corresponding factors in the sum for calculating of the cross section.
203 // Therefore we can just set them equal to 1 in this case. In our opinion it is simplier to eliminate
204 // singularity in such way, because it allows to keep intact original formulas from paper.
205 C_S_plus = C_S_plus==0?1:C_S_plus;
206 C_S_minus = C_S_minus==0?1:C_S_minus;
207
208
209 /* Bkg Contribution */
210 //Auxiliary functions
211 double W_plus = W + M;
212 double W_minus = W - M;
213 double W_plus2 = W_plus*W_plus;
214 double W_minus2 = W_minus*W_minus;
215 double s_minus_M2 = W_plus*W_minus;
216 double u_minus_M2 = m_pi2 - 2*(q_0*p_10 + abs_mom_q*abs_mom_k*CosTheta);
217 double t_minus_mpi2 = -(Q2 + 2*qk);
218 double mpi2_minus_k2 = Q2 + m_pi2;
219 double Q = TMath::Sqrt(Q2);
220
221 // Helicity amplitudes for bkg
222 HelicityBkgAmp<std::complex<double>> Hbkg;
223
224 //Vector_helicity amplituds for bkg
225 // vector cut
226 double FV_cut = 0;
227 //virtual form symmetry_factor to kill the bkg smothly
228 if (W < fBkgVWmin)
229 FV_cut = 1;
230 else if (W >= fBkgVWmin && W < fBkgVWmax)
231 FV_cut = fBkgV0 + W*(fBkgV1 + W*(fBkgV2 + W*fBkgV3));
232
233 fFormFactors.Calculate(interaction);
234 double g_A = -FA0;
235 double FF_pi = fFormFactors.F1V();
236 double FF_1 = 0.5*FF_pi;
237 double FF_2 = 0.5*fFormFactors.xiF2V();
238
239
240
241 // Eq. 38 and Table 5 of ref. 1
242 double V_1 = 0, V_2 = 0, V_3 = 0, V_4 = 0, V_5 = 0;
243
244 // [p,n,pi+,pi0,pi-]
245 switch (spp_channel) {
246
247 //CC
248 case (kSpp_vp_cc_10100) :
249 case (kSpp_vbn_cc_01001):
250 V_1 = g_A*kSqrt2*FV_cut/f_pi*(Mt2*FF_1/u_minus_M2 + FF_2/M);
251 V_2 = g_A*kSqrt2*FV_cut/f_pi*Mt2*FF_1/u_minus_M2/qk;
252 V_3 = g_A*kSqrt2*FV_cut/f_pi*2*FF_2/u_minus_M2;
253 V_4 = -V_3;
254 V_5 = g_A*kSqrt2*FV_cut/f_pi*M*FF_pi/t_minus_mpi2/qk;
255 break;
256
257 case (kSpp_vn_cc_01100) :
258 case (kSpp_vbp_cc_10001):
259 V_1 = g_A*kSqrt2*FV_cut/f_pi*(Mt2*FF_1/s_minus_M2 + FF_2/M);
260 V_2 = g_A*kSqrt2*FV_cut/f_pi*Mt2*FF_1/s_minus_M2/qk;
261 V_3 = -g_A*kSqrt2*FV_cut/f_pi*2*FF_2/s_minus_M2;
262 V_4 = V_3;
263 V_5 = -g_A*kSqrt2*FV_cut/f_pi*M*FF_pi/t_minus_mpi2/qk;
264 break;
265
266 case (kSpp_vn_cc_10010) :
267 case (kSpp_vbp_cc_01010):
268 V_1 = g_A*FV_cut/f_pi*Mt2*FF_1*(1/s_minus_M2 - 1/u_minus_M2);
269 V_2 = g_A*FV_cut/f_pi*Mt2*FF_1/qk*(1/s_minus_M2 - 1/u_minus_M2);
270 V_3 = -g_A*FV_cut/f_pi*2*FF_2*(1/s_minus_M2 + 1/u_minus_M2);
271 V_4 = -g_A*FV_cut/f_pi*2*FF_2*(1/s_minus_M2 - 1/u_minus_M2);;
272 V_5 = -g_A*FV_cut/f_pi*Mt2*FF_pi/t_minus_mpi2/qk;
273 break;
274
275 //NC
276 case (kSpp_vp_nc_10010) :
277 case (kSpp_vbp_nc_10010):
278 case (kSpp_vn_nc_01010) :
279 case (kSpp_vbn_nc_01010):
280 V_1 = g_A*FV_cut*M/f_pi*(FF_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_2/M2);
281 V_2 = g_A*FV_cut*M/f_pi*FF_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
282 V_3 = g_A*FV_cut/2/f_pi*FF_2*(1/u_minus_M2 - 1/s_minus_M2);
283 V_4 = -g_A*FV_cut/2/f_pi*FF_2*(1/u_minus_M2 + 1/s_minus_M2);
284 break;
285
286 case (kSpp_vp_nc_01100) :
287 case (kSpp_vbp_nc_01100):
288 V_1 = -k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2);
289 V_2 = -k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2)/qk;
290 V_3 = -k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 + 1/s_minus_M2);
291 V_4 = k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 - 1/s_minus_M2);
292 V_5 = -k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_pi/t_minus_mpi2/qk;
293 break;
294
295 case (kSpp_vn_nc_10001) :
296 case (kSpp_vbn_nc_10001):
297 V_1 = k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2);
298 V_2 = k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_1*(1/u_minus_M2 - 1/s_minus_M2)/qk;
299 V_3 = k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 + 1/s_minus_M2);
300 V_4 = -k1_Sqrt2*g_A*FV_cut/f_pi*FF_2*(1/u_minus_M2 - 1/s_minus_M2);
301 V_5 = k1_Sqrt2*g_A*FV_cut*Mt2/f_pi*FF_pi/t_minus_mpi2/qk;
302 break;
303 default:
304 // should not be here - meaningless to return anything
305 gAbortingInErr = true;
306 LOG("MKSPPPXSec2020", pFATAL) << "Unknown resonance production channel: " << spp_channel;
307 exit(1);
308
309
310 }
311 double V_25 = (W_plus*W_minus)*V_2 - Q2*V_5;
312
313
314 double O_1_plus = TMath::Sqrt((W_plus2 + Q2)*( W_plus2 - m_pi2 ))/Wt2;
315 // There is no check of positivity under root expression in the original code
316 double O_1_minus = TMath::Sqrt(TMath::Max(0., (W_minus2 + Q2)*(W_minus2 - m_pi2)))/Wt2;
317 double O_2_plus = TMath::Sqrt((W_plus2 + Q2)/(W_plus2 - m_pi2));
318 // There is no check of positivity under root expression in the original code
319 double O_2_minus = W_minus2>m_pi2?TMath::Sqrt((W_minus2 + Q2)/(W_minus2 - m_pi2)):0;
320
321 double K_1_V = W_minus*O_1_plus;
322 double K_2_V = W_plus*O_1_minus;
323 // There is no check of positivity under root expression in the original code
324 double K_3_V = W_minus2>m_pi2?abs_mom_q*abs_mom_q*W_plus*O_2_minus:0;
325 double K_4_V = abs_mom_q*abs_mom_q*W_minus*O_2_plus;
326 double K_5_V = 1/O_2_plus;
327 // the of K_6_V = 1/O_2_minus differ from original code
328 double K_6_V = TMath::Sqrt(TMath::Max(0., (W_minus2 - m_pi2))/(W_minus2 + Q2));
329
330 double F_1 = V_1 + qk*(V_3 - V_4)/W_minus + W_minus*V_4;
331 double F_2 = -V_1 + qk*(V_3 - V_4)/W_plus + W_plus *V_4;
332 double F_3 = (V_3 - V_4) + V_25/W_plus;
333 double F_4 = (V_3 - V_4) - V_25/W_minus;
334 double F_5 = (W_plus2 + Q2)*(V_1 + W_minus*V_4)/Wt2 + (W_plus*q_0 - qk)*(V_3 - V_4) + q_0*V_25 - k_0*qk*V_5 -
335 qk*(W_plus2 + Q2 + 2*W*W_minus)*V_2/Wt2;
336 double F_6 =-(W_minus2 + Q2)*(V_1 - W_plus *V_4)/Wt2 + (W_minus*q_0 - qk)*(V_3 - V_4) - q_0*V_25 + k_0*qk*V_5 +
337 qk*(W_plus2 + Q2 + 2*W*W_minus)*V_2/Wt2;
338
339 double sF_1 = K_1_V*F_1/Mt2;
340 double sF_2 = K_2_V*F_2/Mt2;
341 double sF_3 = K_3_V*F_3/Mt2;
342 double sF_4 = K_4_V*F_4/Mt2;
343 double sF_5 = K_5_V*F_5/Mt2;
344 double sF_6 = K_6_V*F_6/Mt2;
345
346 Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
347 k1_Sqrt2*SinTheta*CosHalfTheta*(sF_3 + sF_4)*std::complex<double>
348 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
349 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
350
351 Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
352 -k1_Sqrt2*SinTheta*SinHalfTheta*(sF_3 - sF_4)*std::complex<double>
353 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
354 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
355
356 Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
357 kSqrt2*(CosHalfTheta*(sF_1 - sF_2) - 0.5*SinTheta*SinHalfTheta*(sF_3 - sF_4))*std::complex<double>
358 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
359 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
360
361 Hbkg(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
362 -kSqrt2*(SinHalfTheta*(sF_1 + sF_2) + 0.5*SinTheta*CosHalfTheta*(sF_3 + sF_4))*std::complex<double>
363 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
364 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
365
366
367 Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
368 -kSqrt2*(SinHalfTheta*(sF_1 + sF_2) + 0.5*SinTheta*CosHalfTheta*(sF_3 + sF_4))*std::complex<double>
369 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
370 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
371
372 Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
373 -kSqrt2*(CosHalfTheta*(sF_1 - sF_2) - 0.5*SinTheta*SinHalfTheta*(sF_3 - sF_4))*std::complex<double>
374 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
375 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
376
377 Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
378 k1_Sqrt2*SinTheta*SinHalfTheta*(sF_3 - sF_4)*std::complex<double>
379 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
380 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
381
382 Hbkg(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
383 k1_Sqrt2*SinTheta*CosHalfTheta*(sF_3 + sF_4)*std::complex<double>
384 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
385 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
386
387
388 if (is_CC || (is_NC && is_nu) )
389 {
390 Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
391 CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 + sF_6)/C_S_minus*std::complex<double>
392 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
393 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
394
395 Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
396 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 - sF_6)/C_S_minus*std::complex<double>
397 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
398 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
399
400 Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
401 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 - sF_6)/C_S_minus*std::complex<double>
402 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
403 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
404
405 Hbkg(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
406 -CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sF_5 + sF_6)/C_S_minus*std::complex<double>
407 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
408 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
409 }
410 if (is_CC || (is_NC && is_nubar) )
411 {
412 Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
413 CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 + sF_6)/C_S_plus*std::complex<double>
414 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
415 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
416
417 Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
418 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 - sF_6)/C_S_plus*std::complex<double>
419 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
420 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
421
422 Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
423 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 - sF_6)/C_S_plus*std::complex<double>
424 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
425 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
426
427 Hbkg(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
428 -CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sF_5 + sF_6)/C_S_plus*std::complex<double>
429 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
430 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
431 }
432
433 if (is_NC)
434 {
435 // isoscalar electromagnetic amplitudes eq.68 of ref. 8
436 fEMFormFactors.Calculate(interaction);
437 double FF_em_1 = 0.5*fEMFormFactors.F1V();
438 double FF_em_2 = 0.5*fEMFormFactors.xiF2V();
439
440 double Vem_1 = 0, Vem_2 = 0, Vem_3 = 0, Vem_4 = 0, Vem_5 = 0;
441
442 // [p,n,pi+,pi0,pi-]
443 switch (spp_channel) {
444 //NC
445 case (kSpp_vp_nc_10010) :
446 case (kSpp_vbp_nc_10010):
447 Vem_1 = -k1_Sqrt2*g_A*FV_cut*M/f_pi*(FF_em_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_em_2/M2);
448 Vem_2 = -k1_Sqrt2*g_A*FV_cut*M/f_pi*FF_em_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
449 Vem_3 = -k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 - 1/s_minus_M2);
450 Vem_4 = k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 + 1/s_minus_M2);
451 break;
452 case (kSpp_vn_nc_01010) :
453 case (kSpp_vbn_nc_01010):
454 Vem_1 = k1_Sqrt2*g_A*FV_cut*M/f_pi*(FF_em_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_em_2/M2);
455 Vem_2 = k1_Sqrt2*g_A*FV_cut*M/f_pi*FF_em_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
456 Vem_3 = k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 - 1/s_minus_M2);
457 Vem_4 =-k1_Sqrt2*g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 + 1/s_minus_M2);
458 break;
459
460 case (kSpp_vp_nc_01100) :
461 case (kSpp_vbp_nc_01100):
462 case (kSpp_vn_nc_10001) :
463 case (kSpp_vbn_nc_10001):
464 Vem_1 = g_A*FV_cut*M/f_pi*(FF_em_1*(1/u_minus_M2 + 1/s_minus_M2) + FF_em_2/M2/2);
465 Vem_2 = g_A*FV_cut*M/f_pi*FF_em_1*(1/u_minus_M2 + 1/s_minus_M2)/qk;
466 Vem_3 = g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 - 1/s_minus_M2);
467 Vem_4 =-g_A*FV_cut/2/f_pi*FF_em_2*(1/u_minus_M2 + 1/s_minus_M2);
468 break;
469 default:
470 // should not be here - meaningless to return anything
471 gAbortingInErr = true;
472 LOG("MKSPPPXSec2020", pFATAL) << "CC channel in NC mode " << spp_channel;
473 exit(1);
474 }
475 double Vem_25 = (W_plus*W_minus)*Vem_2 + Vem_5;
476
477 double Fem_1 = Vem_1 + qk*(Vem_3 - Vem_4)/W_minus + W_minus*Vem_4;
478 double Fem_2 = -Vem_1 + qk*(Vem_3 - Vem_4)/W_plus + W_plus*Vem_4;
479 double Fem_3 = (Vem_3 - Vem_4) + Vem_25/W_plus;
480 double Fem_4 = (Vem_3 - Vem_4) - Vem_25/W_minus;
481 double Fem_5 = (W_plus2 + Q2)*(Vem_1 + W_minus*Vem_4)/Wt2 + (W_plus*q_0 - qk)*(Vem_3 - Vem_4) + q_0*Vem_25 -
482 qk*(W_plus2 + Q2 + 2*W*W_minus)*Vem_2/Wt2;
483 double Fem_6 =-(W_minus2 + Q2)*(Vem_1 - W_plus*Vem_4)/Wt2 + (W_minus*q_0 - qk)*(Vem_3 - Vem_4) - q_0*Vem_25 +
484 qk*(W_plus2 + Q2 + 2*W*W_minus)*Vem_2/Wt2;
485
486 double sFem_1 = K_1_V*Fem_1/Mt2;
487 double sFem_2 = K_2_V*Fem_2/Mt2;
488 double sFem_3 = K_3_V*Fem_3/Mt2;
489 double sFem_4 = K_4_V*Fem_4/Mt2;
490 double sFem_5 = K_5_V*Fem_5/Mt2;
491 double sFem_6 = K_6_V*Fem_6/Mt2;
492
493 HelicityBkgAmp<std::complex<double>> HbkgEM;
494
495 HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
496 k1_Sqrt2*SinTheta*CosHalfTheta*(sFem_3 + sFem_4)*std::complex<double>
497 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
498 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
499
500 HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
501 -k1_Sqrt2*SinTheta*SinHalfTheta*(sFem_3 - sFem_4)*std::complex<double>
502 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
503 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
504
505 HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
506 kSqrt2*(CosHalfTheta*(sFem_1 - sFem_2) - 0.5*SinTheta*SinHalfTheta*(sFem_3 - sFem_4))*std::complex<double>
507 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
508 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
509
510 HbkgEM(Current::VECTOR, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
511 -kSqrt2*(SinHalfTheta*(sFem_1 + sFem_2) + 0.5*SinTheta*CosHalfTheta*(sFem_3 + sFem_4))*std::complex<double>
512 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
513 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
514
515
516 HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
517 -kSqrt2*(SinHalfTheta*(sFem_1 + sFem_2) + 0.5*SinTheta*CosHalfTheta*(sFem_3 + sFem_4))*std::complex<double>
518 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
519 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
520
521 HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
522 -kSqrt2*(CosHalfTheta*(sFem_1 - sFem_2) - 0.5*SinTheta*SinHalfTheta*(sFem_3 - sFem_4))*std::complex<double>
523 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
524 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
525
526 HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
527 k1_Sqrt2*SinTheta*SinHalfTheta*(sFem_3 - sFem_4)*std::complex<double>
528 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
529 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
530
531 HbkgEM(Current::VECTOR, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
532 k1_Sqrt2*SinTheta*CosHalfTheta*(sFem_3 + sFem_4)*std::complex<double>
533 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
534 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
535
536
537 if (is_nu)
538 {
539 HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
540 CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 + sFem_6)/C_S_minus*std::complex<double>
541 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
542 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
543
544 HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
545 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 - sFem_6)/C_S_minus*std::complex<double>
546 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
547 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
548
549 HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
550 -SinHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 - sFem_6)/C_S_minus*std::complex<double>
551 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
552 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
553
554 HbkgEM(Current::VECTOR, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
555 -CosHalfTheta*(k_0*eps_z_L - abs_mom_k*eps_zero_L)*(sFem_5 + sFem_6)/C_S_minus*std::complex<double>
556 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
557 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
558 }
559 else
560 {
561 HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
562 CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 + sFem_6)/C_S_plus*std::complex<double>
563 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
564 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
565
566 HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
567 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 - sFem_6)/C_S_plus*std::complex<double>
568 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
569 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
570
571 HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
572 -SinHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 - sFem_6)/C_S_plus*std::complex<double>
573 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
574 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
575
576 HbkgEM(Current::VECTOR, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
577 -CosHalfTheta*(k_0*eps_z_R - abs_mom_k*eps_zero_R)*(sFem_5 + sFem_6)/C_S_plus*std::complex<double>
578 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
579 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
580 }
581
582 Hbkg = (1 - 2*fSin2Wein)*Hbkg - 4*fSin2Wein*HbkgEM;
583 }
584
585 //Axial helicity amp for bkg
586 //Axial cut is same as FV_cut in the latest version
587 double FA_cut = FV_cut;
588
589 // FF_A here is equal to -g_A*FF_A of original code
590 double FF_A = fFormFactors.FA();
591 double F_rho = Frho0/(1 - (t_minus_mpi2 + m_pi2 )/fRho770Mass/fRho770Mass);
592
593 // Eq. 38 and Table 5 of ref. 1
594 double A_1 = 0, A_3 = 0, A_4 = 0, A_7 = 0, A_8 = 0;
595
596 // [p,n,pi+,pi0,pi-]
597 switch (spp_channel) {
598
599 //CC
600 case (kSpp_vp_cc_10100) :
601 case (kSpp_vbn_cc_01001):
602 A_1 = -kSqrt2*FA_cut*M*FF_A/f_pi/u_minus_M2;
603 A_3 = -A_1;
604 A_4 = k1_Sqrt2*FA_cut/f_pi*(FF_A + F_rho)/M;
605 A_7 = kSqrt2*FA_cut/f_pi*M*FF_A/mpi2_minus_k2;
606 A_8 = -k1_Sqrt2*FA_cut/f_pi*(FF_A*(1 + 4*M2/u_minus_M2) + F_rho)/mpi2_minus_k2;
607 break;
608
609 case (kSpp_vn_cc_01100) :
610 case (kSpp_vbp_cc_10001):
611 A_1 = kSqrt2*FA_cut/f_pi*M*FF_A/s_minus_M2;
612 A_3 = A_1;
613 A_4 = -k1_Sqrt2*FA_cut/f_pi/M*(FF_A + F_rho);
614 A_7 = kSqrt2*FA_cut/f_pi*M*FF_A/mpi2_minus_k2;
615 A_8 = k1_Sqrt2/f_pi*FA_cut*(FF_A/mpi2_minus_k2*(1 + 4*M2/s_minus_M2) + F_rho/mpi2_minus_k2);
616 break;
617
618 case (kSpp_vn_cc_10010) :
619 case (kSpp_vbp_cc_01010):
620 A_1 = FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2);
621 A_3 = -FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2) ;
622 A_4 = -FA_cut/f_pi/M*(FF_A + F_rho);
623 A_8 = FA_cut/f_pi*(F_rho + FF_A*(1 + 2*M2*(1/u_minus_M2 + 1/s_minus_M2)))/mpi2_minus_k2;
624 break;
625
626 //NC
627 case (kSpp_vp_nc_10010) :
628 case (kSpp_vbp_nc_10010):
629 case (kSpp_vn_nc_01010) :
630 case (kSpp_vbn_nc_01010):
631 A_1 = -FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2)/2;
632 A_3 = FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2)/2;
633 A_7 = FA_cut/f_pi*M*FF_A/mpi2_minus_k2;
634 break;
635
636 case (kSpp_vp_nc_01100) :
637 case (kSpp_vbp_nc_01100):
638 A_1 = k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2);
639 A_3 = -k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2);
640 A_4 = -k1_Sqrt2*FA_cut/f_pi/M*(FF_A + F_rho);
641 A_8 = k1_Sqrt2*FA_cut/f_pi/mpi2_minus_k2*(F_rho + FF_A*(1 + 2*M2*(1/u_minus_M2 + 1/s_minus_M2)));
642 break;
643
644 case (kSpp_vn_nc_10001) :
645 case (kSpp_vbn_nc_10001):
646 A_1 = -k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 + 1/s_minus_M2);
647 A_3 = k1_Sqrt2*FA_cut/f_pi*M*FF_A*(1/u_minus_M2 - 1/s_minus_M2);
648 A_4 = k1_Sqrt2*FA_cut/f_pi/M*(FF_A + F_rho);
649 A_8 = k1_Sqrt2*FA_cut/f_pi/mpi2_minus_k2*(FF_A*(1 + 2*M2*(1/u_minus_M2 + 1/s_minus_M2)) + F_rho);
650 break;
651 default:
652 // should not be here - meaningless to return anything
653 gAbortingInErr = true;
654 LOG("MKSPPPXSec2020", pFATAL) << "Unknown resonance production channel: " << spp_channel;
655 exit(1);
656 }
657
658 double K_1_A = abs_mom_q*O_2_plus;
659 // There is no check of positivity under root expression in the original code
660 double K_2_A = W_minus2>m_pi2?abs_mom_q*O_2_minus:1;
661 double K_3_A = abs_mom_q*O_1_minus;
662 double K_4_A = abs_mom_q* O_1_plus;
663 double K_5_A = O_1_minus;
664 double K_6_A = O_1_plus;
665 double abs_mom_k2 = abs_mom_k*abs_mom_k;
666 double Delta = k_0*(q_0*k_0 - qk)/abs_mom_k2;
667
668 double G_1 = W_plus* A_1 - M*A_4;
669 double G_2 = -W_minus*A_1 - M*A_4;
670 double G_3 = A_1 - A_3;
671 double G_4 = -G_3;
672 double G_5 = (Delta + (W_plus2 - m_pi2)/Wt2 + Wt2*k_0*W_plus/(W_minus2 + Q2))*A_1 + (q_0 - Delta)*A_3 - M*W_minus*A_4/(p_10 - M);
673 double G_6 =-(Delta + (W_minus2 - m_pi2)/Wt2 + Wt2*k_0*W_minus/(W_plus2 + Q2))*A_1 - (q_0 - Delta)*A_3 - M*W_plus*A_4 /(p_10 + M);
674 double G_7= (W_plus2 - m_pi2)*A_1/Wt2 + q_0*A_3 - M*A_4 + k_0*A_7 + W_plus *k_0*A_8;
675 double G_8 = -(W_minus2 - m_pi2)*A_1/Wt2 - q_0*A_3 - M*A_4 - k_0*A_7 + W_minus*k_0*A_8;
676
677 //Isobaric amplitud (Scribal G)
678 double sG_1 = K_1_A*G_1/Mt2;
679 double sG_2 = K_2_A*G_2/Mt2;
680 double sG_3 = K_3_A*G_3/Mt2;
681 double sG_4 = K_4_A*G_4/Mt2;
682 double sG_5 = K_5_A*G_5/Mt2;
683 double sG_6 = K_6_A*G_6/Mt2;
684 double sG_7 = K_5_A*G_7/Mt2;
685 double sG_8 = K_6_A*G_8/Mt2;
686
687 // scalar part eq. 71 of ref.8
688 if (is_NC)
689 {
690 double g_s = 0.15;
691 if (spp_channel == kSpp_vp_nc_01100 || spp_channel == kSpp_vbp_nc_01100 || spp_channel == kSpp_vn_nc_10001 || spp_channel == kSpp_vbn_nc_10001)
692 g_s *= kSqrt2;
693
694 double Gs_A = g_s/g_A*FF_A;
695 double As_1 = k1_Sqrt2*g_A*FA_cut/f_pi*M*Gs_A*(1/u_minus_M2 - 1/s_minus_M2);
696 double As_3 =-k1_Sqrt2*g_A*FA_cut/f_pi*M*Gs_A*(1/u_minus_M2 + 1/s_minus_M2);
697 double As_4 = 0;
698 double As_7 = 0;
699 double As_8 = 0;
700 if (spp_channel == kSpp_vn_nc_01010 || spp_channel == kSpp_vbn_nc_01010)
701 {
702 As_1 *= -1;
703 As_3 *= -1;
704 }
705
706 double Gs_1 = W_plus* As_1 - M*As_4;
707 double Gs_2 = -W_minus*As_1 - M*As_4;
708 double Gs_3 = As_1 - As_3;
709 double Gs_4 = -Gs_3;
710 double Gs_5 = (Delta + (W_plus2 - m_pi2)/Wt2 + 2*W*k_0*W_plus /(W_minus2 + Q2))*As_1 + (q_0 - Delta)*As_3 - M*W_minus*As_4/(p_10-M);
711 double Gs_6 =-(Delta + (W_minus2 - m_pi2)/Wt2 + 2*W*k_0*W_minus/(W_plus2 + Q2)) *As_1 - (q_0 - Delta)*As_3 - M*W_plus*As_4/(p_10 + M);
712 double Gs_7 = (W_plus2 - m_pi2)*As_1/Wt2 + q_0*As_3 - M*As_4 + k_0*As_7 + W_plus* k_0*As_8;
713 double Gs_8 =-(W_minus2 - m_pi2)*As_1/Wt2 - q_0*As_3 - M*As_4 - k_0*As_7 + W_minus* k_0*As_8;
714
715 sG_1 -= K_1_A*Gs_1/Mt2;
716 sG_2 -= K_2_A*Gs_2/Mt2;
717 sG_3 -= K_3_A*Gs_3/Mt2;
718 sG_4 -= K_4_A*Gs_4/Mt2;
719 sG_5 -= K_5_A*Gs_5/Mt2;
720 sG_6 -= K_6_A*Gs_6/Mt2;
721 sG_7 -= K_5_A*Gs_7/Mt2;
722 sG_8 -= K_6_A*Gs_8/Mt2;
723 }
724
725 // helicity amplitudes
726 Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
727 k1_Sqrt2*SinTheta*CosHalfTheta*(sG_3 + sG_4)*std::complex<double>
728 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
729 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
730
731 Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
732 k1_Sqrt2*SinTheta*SinHalfTheta*(sG_3 - sG_4)*std::complex<double>
733 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
734 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
735
736 Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
737 kSqrt2*(CosHalfTheta*(sG_1 - sG_2) - 0.5*SinTheta*SinHalfTheta*(sG_3 - sG_4))*std::complex<double>
738 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
739 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
740
741 Hbkg(Current::AXIAL, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
742 kSqrt2*(SinHalfTheta*(sG_1 + sG_2) + 0.5*SinTheta*CosHalfTheta*(sG_3 + sG_4))*std::complex<double>
743 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
744 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
745
746
747 Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
748 -kSqrt2*(SinHalfTheta*(sG_1 + sG_2) + 0.5*SinTheta*CosHalfTheta*(sG_3 + sG_4))*std::complex<double>
749 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
750 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
751
752 Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
753 kSqrt2*(CosHalfTheta*(sG_1 - sG_2) - 0.5*SinTheta*SinHalfTheta*(sG_3 - sG_4))*std::complex<double>
754 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
755 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
756
757 Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
758 k1_Sqrt2*SinTheta*SinHalfTheta*(sG_3 - sG_4)*std::complex<double>
759 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
760 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
761
762 Hbkg(Current::AXIAL, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
763 -k1_Sqrt2*SinTheta*CosHalfTheta*(sG_3 + sG_4)*std::complex<double>
764 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
765 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
766
767 if (is_CC || (is_NC && is_nu) )
768 {
769 Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
770 CosHalfTheta*(abs_mom_k*eps_z_L*(sG_5 + sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 + sG_8))/k_0/C_S_minus*std::complex<double>
771 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
772 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
773
774 Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
775 SinHalfTheta*(abs_mom_k*eps_z_L*(sG_5 - sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 - sG_8))/k_0/C_S_minus*std::complex<double>
776 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
777 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
778
779 Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
780 -SinHalfTheta*(abs_mom_k*eps_z_L*(sG_5 - sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 - sG_8))/k_0/C_S_minus*std::complex<double>
781 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
782 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
783
784 Hbkg(Current::AXIAL, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
785 CosHalfTheta*(abs_mom_k*eps_z_L*(sG_5 + sG_6) + (k_0*eps_zero_L - abs_mom_k*eps_z_L)*(sG_7 + sG_8))/k_0/C_S_minus*std::complex<double>
786 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
787 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
788 }
789 if (is_CC || (is_NC && is_nubar) )
790 {
791 Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
792 CosHalfTheta*(abs_mom_k*eps_z_R*(sG_5 + sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 + sG_8))/k_0/C_S_plus*std::complex<double>
793 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
794 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
795
796 Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
797 SinHalfTheta*(abs_mom_k*eps_z_R*(sG_5 - sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 - sG_8))/k_0/C_S_plus*std::complex<double>
798 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
799 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
800
801 Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
802 -SinHalfTheta*(abs_mom_k*eps_z_R*(sG_5 - sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 - sG_8))/k_0/C_S_plus*std::complex<double>
803 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
804 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
805
806 Hbkg(Current::AXIAL, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
807 CosHalfTheta*(abs_mom_k*eps_z_R*(sG_5 + sG_6) + (k_0*eps_zero_R - abs_mom_k*eps_z_R)*(sG_7 + sG_8))/k_0/C_S_plus*std::complex<double>
808 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
809 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
810 }
811
812 /* Resonace contribution */
813 HelicityAmpVminusARes<std::complex<double>> Hres;
814 //f_BW_A is same as f_BW_V in the latest version, so rename f_BW_V = f_BW_A -> f_BW
815 std::complex<double> kappa_f_BW;
816
817 for (auto res : fResList)
818 {
819
820 if (utils::res::Isospin(res) == 1 && SppChannel::FinStateIsospin(spp_channel) == 3) // skip resonances with I=1/2 if isospin of final state is 3/2
821 continue;
822
823 int NR = utils::res::ResonanceIndex (res);
824 int LR = utils::res::OrbitalAngularMom (res);
825 int JR = (utils::res::AngularMom (res) + 1)/2;
826 double MR = utils::res::Mass (res);
827 double WR = utils::res::Width (res);
828 double Cjsgn_plus = utils::res::Cjsgn_plus (res);
829 double Dsgn = utils::res::Dsgn (res);
830 double BR = SppChannel::BranchingRatio (res);
831
832
833 double d = W_plus2 + Q2;
834 double sq2omg = TMath::Sqrt(2/fOmega);
835 double nomg = NR*fOmega;
836
837 //Graczyk and Sobczyk vector form-factors
838 double CV_factor = 1/(1 + Q2/fMv2/4);
839 // Eq. 29 of ref. 6
840 double CV3 = fCv3*CV_factor/(1 + Q2/fMv2)/(1 + Q2/fMv2);
841 // Eq. 30 of ref. 6
842 double CV4 = -1. * fCv4 / fCv3 * CV3;
843 // Eq. 31 of ref. 6
844 double CV5 = fCv51*CV_factor/(1 + Q2/fMv2/fCv52)/(1 + Q2/fMv2/fCv52);
845
846 // Eq. 38 of ref. 6
847 double GV3 = 0.5*k1_Sqrt3*(CV4*(W2 - M2 - Q2)/2/M2 + CV5*(W2 - M2 + Q2)/2/M2 + CV3*W_plus/M);
848 // Eq. 39 of ref. 6
849 double GV1 = -0.5*k1_Sqrt3*(CV4*(W2 - M2 - Q2)/2/M2 + CV5*(W2 - M2 + Q2)/2/M2 - CV3*(W_plus*M + Q2)/W/M);
850 // Eq. 36 of ref. 6, which is implied to use for EM-production
851 double GV = 0.5*TMath::Sqrt(1 + Q2/W_plus2)/TMath::Power(1 + Q2/4/M2, 0.5*NR)*TMath::Sqrt(3*GV3*GV3 + GV1*GV1);
852 // Eq. 37 of ref. 6, which is implied to use for neutrino-production
853 // double GV = 0.5*TMath::Sqrt(1 + Q2/W_plus2)/TMath::Power(1 + Q2/4/M2, NR)*TMath::Sqrt(3*GV3*GV3 + GV1*GV1);
854
855
856 //Graczyk and Sobczyk axial form-symmetry_factor
857 // Eq. 52 of ref. 6
858 double CA5 = fCA50/(1 + Q2/fMa2)/(1 + Q2/fMa2);
859
860 // The form is the same like in Eq. 54 of ref. 6, but differ from it by index, which in ref. 6 is equal to NR.
861 double GA = 0.5*kSqrt3*TMath::Sqrt(1 + Q2/W_plus2)*(1 - (W2 - Q2 -M2)/8/M2)*CA5/TMath::Power(1+ Q2/4/M2, 0.5*NR);
862
863 double qMR_0 = (MR*MR - M2 + m_pi2)/(2*MR);
864 double abs_mom_qMR = TMath::Sqrt( qMR_0*qMR_0 - m_pi2);
865 double Gamma = WR*TMath::Power((abs_mom_q/abs_mom_qMR), 2*LR + 1);
866
867 // denominator of Breit-Wigner function
868 std::complex<double> denom(W - MR, Gamma/2);
869 // Breit-Wigner amplitude multiplied by kappa*sqrt(BR) to avoid singularity at abs_mom_q=0, where BR = chi_E (see eq. 25 and 27 of ref. 1)
870 // double kappa = kPi*W*TMath::Sqrt(2/JR/abs_mom_q)/M;
871 // f_BW = TMath::Sqrt(BR*Gamma/2/kPi)/denom;
872 kappa_f_BW = W*TMath::Sqrt(kPi*BR*WR/JR/abs_mom_qMR)*TMath::Power((abs_mom_q/abs_mom_qMR), LR)/denom/M;
873
874 fFKR.Lamda = sq2omg*abs_mom_k;
875 fFKR.Tv = GV/3/W/sq2omg;
876 fFKR.Ta = 2./3/sq2omg*abs_mom_k*GA/d;
877 fFKR.Rv = kSqrt2*abs_mom_k*W_plus*GV/d;
878 fFKR.Ra = kSqrt2/6*(W_plus + 2*nomg*W/d)*GA/W;
879 fFKR.R = fFKR.Rv;
880 fFKR.T = fFKR.Tv;
881 fFKR.Rplus = - (fFKR.Rv + fFKR.Ra);
882 fFKR.Rminus = - (fFKR.Rv - fFKR.Ra);
883 fFKR.Tplus = - (fFKR.Tv + fFKR.Ta);
884 fFKR.Tminus = - (fFKR.Tv - fFKR.Ta);
885
886 double a_aux = 1 + ((W2 + Q2 + M2)/(Mt2*W));
887 // if Q2=Q=sqrt(Q2)=0 then the singularity appears in k_sqrtQ2
888 // which is eliminated when in Hres at BosonPolarization=PLUS0 or BosonPolarization=MINUS0
889 // when k_sqrtQ2 is multiplied by C, B and S. Therefore in this case we put Q equal to 1.
890 // In our opinion it is simplier to eliminate singularity in such way, because it allows to
891 // keep intact original formulas from paper.
892 if (Q == 0) Q = 1;
893
894 double C_plus = Q/C_S_plus*((eps_zero_R*abs_mom_k - eps_z_R*k_0)*(1./3 + k_0/a_aux/M) +
895 (Wt2/3 - Q2/a_aux/M + nomg/a_aux/M/3)*(eps_z_R + (eps_zero_R*k_0 - eps_z_R*abs_mom_k)*abs_mom_k/mpi2_minus_k2))*GA/Wt2/abs_mom_k;
896 double C_minus = Q/C_S_minus*((eps_zero_L*abs_mom_k - eps_z_L*k_0)*(1./3 + k_0/a_aux/M) +
897 (Wt2/3 - Q2/a_aux/M + nomg/a_aux/M/3)*(eps_z_L + (eps_zero_L*k_0 - eps_z_L*abs_mom_k)*abs_mom_k/mpi2_minus_k2))*GA/Wt2/abs_mom_k;
898
899 double B_plus = Q/C_S_plus *(eps_zero_R + eps_z_R*abs_mom_k/a_aux/M +
900 (eps_zero_R*k_0 - eps_z_R*abs_mom_k)*(k_0 + abs_mom_k*abs_mom_k/M/a_aux)/mpi2_minus_k2)*GA/W/3/sq2omg/abs_mom_k;
901 double B_minus = Q/C_S_minus*(eps_zero_L + eps_z_L*abs_mom_k/a_aux/M +
902 (eps_zero_L*k_0 - eps_z_L*abs_mom_k)*(k_0 + abs_mom_k*abs_mom_k/M/a_aux)/mpi2_minus_k2)*GA/W/3/sq2omg/abs_mom_k;
903
904 double S_plus = Q/C_S_plus *(eps_z_R*k_0 - eps_zero_R*abs_mom_k)*(1 + Q2/M2 - 3*W/M)*GV/abs_mom_k_L2/6;
905 double S_minus = Q/C_S_minus*(eps_z_L*k_0 - eps_zero_L*abs_mom_k)*(1 + Q2/M2 - 3*W/M)*GV/abs_mom_k_L2/6;
906 double k_sqrtQ2 = abs_mom_k/Q;
907
908 const RSHelicityAmplModelI * hamplmod = 0;
909
910 if (is_CC)
911 hamplmod = fHAmplModelCC;
912 else if(is_NC)
913 {
914 if (is_p)
915 hamplmod = fHAmplModelNCp;
916 else
917 hamplmod = fHAmplModelNCn;
918 }
919
920 fFKR.S = S_minus;
921 fFKR.B = B_minus;
922 fFKR.C = C_minus;
923 const RSHelicityAmpl & hampl = hamplmod->Compute(res, fFKR);
924 double fp3 = hampl.AmpPlus3();
925 double fp1 = hampl.AmpPlus1();
926 double fm3 = hampl.AmpMinus3();
927 double fm1 = hampl.AmpMinus1();
928 double fm0m = 0, fm0p = 0, fp0m = 0, fp0p = 0;
929 if (is_CC || (is_NC && is_nu) )
930 {
931 fm0m = hampl.Amp0Minus();
932 fm0p = hampl.Amp0Plus();
933 }
934 if (is_CC || (is_NC && is_nubar) )
935 {
936 fFKR.S = S_plus;
937 fFKR.B = B_plus;
938 fFKR.C = C_plus;
939 const RSHelicityAmpl & hampl_plus = hamplmod->Compute(res, fFKR);
940 fp0m = hampl_plus.Amp0Minus();
941 fp0p = hampl_plus.Amp0Plus();
942 }
943
944 double JRtSqrt2 = kSqrt2*JR;
945
946 // The signs of Hres are not correct. Now we consider these signs are exactly equal to ones in the latest version
947 // of code provided by Minoo Kabirnezhad, however pay attention to Cjsgn_plus
948 // which differ from what stated in refs. 1 and 2.
949 // More details about other problems can be found in Ref. 12, see section "Ambiguity in calculation of signs of the amplitudes"
950 // (most important conclusion in subsection "Paradox and probable explanation").
951
952 // Helicity amplitudes V-A, eq. 23 - 25 and Table 3 of ref. 1
953 // Cjsgn_plus are C^j_{lS2z} for S2z=1/2 and given in Table 9 of ref. 4
954 // Cjsgn_minus are C^j_{lS2z} for S2z=-1/2 and all equal to 1 (see Table 7 of ref. 4)
955
956 // The sign of the following amplitude is opposite to one from original code, because of it will be change
957 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
958 Hres(res, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
959 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp3*std::complex<double>
960 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
961 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
962
963 Hres(res, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
964 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp3*std::complex<double>
965 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
966 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
967
968 Hres(res, BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
969 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp1*std::complex<double>
970 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
971 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
972
973
974 // The sign of the following amplitude is opposite to one from original code, because of it will be change
975 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
976 Hres(res, BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
977 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fp1*std::complex<double>
978 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
979 TMath::Sin(Phi*PhaseFactor(BosonPolarization::LEFT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
980
981
982
983 Hres(res, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
984 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm1*std::complex<double>
985 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
986 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
987
988 Hres(res, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
989 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm1*std::complex<double>
990 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
991 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
992
993 Hres(res, BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
994 -JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm3*std::complex<double>
995 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
996 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
997
998 Hres(res, BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
999 JRtSqrt2*Dsgn*Cjsgn_plus*kappa_f_BW*fm3*std::complex<double>
1000 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
1001 TMath::Sin(Phi*PhaseFactor(BosonPolarization::RIGHT, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
1002
1003
1004
1005 Hres(res, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
1006 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0m*std::complex<double>
1007 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
1008 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
1009
1010 // The sign of the following amplitude is opposite to one from original code, because of it will be change
1011 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
1012 Hres(res, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
1013 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0m*std::complex<double>
1014 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
1015 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
1016
1017 Hres(res, BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
1018 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0p*std::complex<double>
1019 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
1020 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
1021
1022 Hres(res, BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
1023 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fm0p*std::complex<double>
1024 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
1025 TMath::Sin(Phi*PhaseFactor(BosonPolarization::MINUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
1026
1027
1028
1029 Hres(res, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS) =
1030 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0m*std::complex<double>
1031 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)),
1032 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::PLUS)));
1033
1034 // The sign of the following amplitude is opposite to one from original code, because of it will be change
1035 // when it will be multiplied by Wigner functions d^j_{\lambda\mu}
1036 Hres(res, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS) =
1037 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0m*std::complex<double>
1038 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)),
1039 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::PLUS)));
1040
1041 Hres(res, BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS) =
1042 k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0p*std::complex<double>
1043 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)),
1044 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::PLUS, NucleonPolarization::MINUS)));
1045
1046 Hres(res, BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS) =
1047 -k_sqrtQ2*JRtSqrt2*Dsgn*kappa_f_BW*fp0p*std::complex<double>
1048 (TMath::Cos(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)),
1049 TMath::Sin(Phi*PhaseFactor(BosonPolarization::PLUS0, NucleonPolarization::MINUS, NucleonPolarization::MINUS)));
1050 } //end resonances loop
1051
1052 // Sum of all helicities
1053 SumHelicityAmpVminusARes<std::complex<double>> sum3;
1054 SumHelicityAmpVminusARes<std::complex<double>> sum1;
1055
1056 double pt_1 = 1;
1057 double pt_2 = 3*CosTheta;
1058 double pt_3 = 15./2*CosTheta*CosTheta - 3./2;
1059 double pt_4 = 35./2*CosTheta*CosTheta*CosTheta - 15./2*CosTheta;
1060
1061 // Wigner functions in terms of pion polar angle (eq. B4 of ref. 1)
1062 double d[4][2][2];
1063 d[0][0][1] = CosHalfTheta;
1064 d[0][0][0] =-SinHalfTheta;
1065 d[0][1][1] = 0;
1066 d[0][1][0] = 0;
1067
1068 d[1][0][1] = CosHalfTheta*(pt_2 - pt_1)/2;
1069 d[1][0][0] =-SinHalfTheta*(pt_2 + pt_1)/2;
1070 d[1][1][1] =-SinHalfTheta*(k1_Sqrt3*pt_2 + kSqrt3*pt_1)/2;
1071 d[1][1][0] = CosHalfTheta*(-k1_Sqrt3*pt_2 + kSqrt3*pt_1)/2;
1072
1073 d[2][0][1] = CosHalfTheta*(pt_3 - pt_2)/3;
1074 d[2][0][0] =-SinHalfTheta*(pt_3 + pt_2)/3;
1075 d[2][1][1] =-SinHalfTheta*(k1_Sqrt2*pt_3 + kSqrt2*pt_2)/3;
1076 d[2][1][0] = CosHalfTheta*(-k1_Sqrt2*pt_3 + kSqrt2*pt_2)/3;
1077
1078 d[3][0][1] = CosHalfTheta*(pt_4 - pt_3)/4;
1079 d[3][0][0] =-SinHalfTheta*(pt_4 + pt_3)/4;
1080 d[3][1][1] =-SinHalfTheta*(kSqrt3_5*pt_4 + kSqrt5_3*pt_3)/4;
1081 d[3][1][0] = CosHalfTheta*(-kSqrt3_5*pt_4 + kSqrt5_3*pt_3)/4;
1082
1083
1084 for (BosonPolarization lk : BosonPolarizationIterator() )
1085 {
1086 int lambda_k = Lambda(lk);
1087 for (NucleonPolarization l2 : NucleonPolarizationIterator() )
1088 {
1089 int lambda_2 = Lambda(l2);
1090 for (NucleonPolarization l1 : NucleonPolarizationIterator() )
1091 {
1092 int lambda_1 = Lambda(l1);
1093 int lambda = lambda_k - lambda_1;
1094 int mu = -lambda_2;
1095 // Symmetry properties of Wigner d-functions, see eq. A1 from ref. 7
1096 int symmetry_factor = 1;
1097 int itemp;
1098 if (lambda < 0)
1099 {
1100 if (TMath::Abs(lambda)>TMath::Abs(mu)) //swap 1 + swap 2
1101 {
1102 symmetry_factor = TMath::Power(-1, (lambda - mu)/2);
1103 lambda*=-1;
1104 mu*=-1;
1105 }
1106 else
1107 {
1108 if (mu<0) //swap 1
1109 {
1110 itemp = lambda;
1111 lambda = -mu;
1112 mu = -itemp;
1113 }
1114 else //swap 2
1115 {
1116 symmetry_factor = TMath::Power(-1, (lambda - mu)/2);
1117 itemp = lambda;
1118 lambda = mu;
1119 mu = itemp;
1120 }
1121 }
1122 }
1123
1124 for (auto res : fResList)
1125 {
1126 // Get baryon resonance parameters
1127 int IR = utils::res::Isospin (res);
1128 int JR = (utils::res::AngularMom (res) + 1)/2;
1129 if (SppChannel::FinStateIsospin(spp_channel) == 3 && IR == 1) // skip resonances with I=1/2 if isospin of final state is 3/2
1130 continue;
1131 // Eq. 24 of ref. 1
1132 if (IR == 3)
1133 sum3(lk, l2, l1) += symmetry_factor*d[JR - 1][(lambda - 1)/2][(mu + 1)/2]*Hres(res, lk, l2, l1);
1134
1135
1136 if (IR == 1)
1137 sum1(lk, l2, l1) += symmetry_factor*d[JR - 1][(lambda - 1)/2][(mu + 1)/2]*Hres(res, lk, l2, l1);
1138 }
1139 }
1140 }
1141 }
1142
1143 // Isospin Clebsch–Gordan coefficients to sum amplitudes for I=1/2 and I=3/2, see eq.25 and Table 2 of ref. 1
1144 double C1 = SppChannel::Isospin1Coefficients(spp_channel);
1145 double C3 = SppChannel::Isospin3Coefficients(spp_channel);
1146
1147
1148 double g = fFermiConstant ;
1149 if(is_CC) g *= fVud;
1150
1151 double Lcoeff= abs_mom_k_L/2/kSqrt2/E;
1152 // Eq. 3.59 of ref. 2, which is multiplied by Q to avoid singularity at Q2=0
1153 double xsec0 = TMath::Power(g/8/kPi/kPi, 2)*abs_mom_q*Lcoeff*Lcoeff/abs_mom_k_L2/2;
1154 // We divide xsec0 by Q2 due to redifinition of Lcoeff
1155
1156
1157 if (kpsdim == 4)
1158 {
1159 for (NucleonPolarization l2 : NucleonPolarizationIterator() )
1160 {
1161 for (NucleonPolarization l1 : NucleonPolarizationIterator() )
1162 {
1163 // Eq. 18 ref. 1
1164 xsec +=
1165 TMath::Power(
1166 std::abs(
1167 C_L_minus*(
1168 Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1) +
1169 C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)
1170 ) +
1171 C_R_minus*(
1172 Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1) +
1173 C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)
1174 ) +
1175 C_S_minus*(
1176 Hbkg(Current::VECTOR, BosonPolarization::MINUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::MINUS0, l2, l1) +
1177 C1*sum1( BosonPolarization::MINUS0, l2, l1) + C3*sum3( BosonPolarization::MINUS0, l2, l1)
1178 )
1179 )
1180 , 2)
1181 +
1182 TMath::Power(
1183 std::abs(
1184 C_L_plus*(
1185 Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1) +
1186 C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)
1187 ) +
1188 C_R_plus*(
1189 Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1) +
1190 C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)
1191 ) +
1192 C_S_plus*(
1193 Hbkg(Current::VECTOR, BosonPolarization::PLUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::PLUS0, l2, l1) +
1194 C1*sum1( BosonPolarization::PLUS0, l2, l1) + C3*sum3( BosonPolarization::PLUS0, l2, l1)
1195 )
1196 )
1197 , 2);
1198 }
1199 }
1200 }
1201 else if (kpsdim == 3)
1202 {
1203 for (NucleonPolarization l2 : NucleonPolarizationIterator() )
1204 {
1205 for (NucleonPolarization l1 : NucleonPolarizationIterator() )
1206 {
1207 xsec +=
1208 (C_L_minus*C_L_minus + C_L_plus*C_L_plus)*(
1209 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1)).real(), 2) +
1210 TMath::Power(std::abs(C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)), 2) +
1211 2*( Hbkg(Current::VECTOR, BosonPolarization::LEFT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::LEFT, l2, l1)).real()*
1212 (C1*sum1( BosonPolarization::LEFT, l2, l1) + C3*sum3( BosonPolarization::LEFT, l2, l1)).real()
1213 )+
1214 (C_R_minus*C_R_minus + C_R_plus*C_R_plus)*(
1215 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1)).real(), 2) +
1216 TMath::Power(std::abs(C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)), 2) +
1217 2*( Hbkg(Current::VECTOR, BosonPolarization::RIGHT, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::RIGHT, l2, l1)).real()*
1218 (C1*sum1( BosonPolarization::RIGHT, l2, l1) + C3*sum3( BosonPolarization::RIGHT, l2, l1)).real()
1219 )+
1220 (C_S_minus*C_S_minus)*(
1221 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::MINUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::MINUS0, l2, l1)).real(), 2) +
1222 TMath::Power(std::abs(C1*sum1( BosonPolarization::MINUS0, l2, l1) + C3*sum3( BosonPolarization::MINUS0, l2, l1)), 2) +
1223 2*( Hbkg(Current::VECTOR, BosonPolarization::MINUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::MINUS0, l2, l1)).real()*
1224 (C1*sum1( BosonPolarization::MINUS0, l2, l1) + C3*sum3( BosonPolarization::MINUS0, l2, l1)).real()
1225 )+
1226 (C_S_plus*C_S_plus)*(
1227 TMath::Power((Hbkg(Current::VECTOR, BosonPolarization::PLUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::PLUS0, l2, l1)).real(), 2) +
1228 TMath::Power(std::abs(C1*sum1( BosonPolarization::PLUS0, l2, l1) + C3*sum3( BosonPolarization::PLUS0, l2, l1)), 2) +
1229 2*( Hbkg(Current::VECTOR, BosonPolarization::PLUS0, l2, l1) - Hbkg(Current::AXIAL, BosonPolarization::PLUS0, l2, l1)).real()*
1230 (C1*sum1( BosonPolarization::PLUS0, l2, l1) + C3*sum3( BosonPolarization::PLUS0, l2, l1)).real()
1231 );
1232 }
1233 }
1234 xsec*= 2*kPi;
1235 }
1236
1237 xsec*=xsec0;
1238
1239
1240 // The algorithm computes d^4xsec/dWdQ2dCosTheta_pidPhi_pi or d^3xsec/dWdQ2dCosTheta_pi
1241 // Check whether variable tranformation is needed
1242 if ( kps != kPSWQ2ctpphipfE && kps != kPSWQ2ctpfE )
1243 {
1244 double J = 1;
1245 if (kpsdim == 3)
1246 J = utils::kinematics::Jacobian(interaction, kPSWQ2ctpfE, kps);
1247 else if (kpsdim == 4)
1248 J = utils::kinematics::Jacobian(interaction, kPSWQ2ctpphipfE, kps);
1249 xsec *= J;
1250 }
1251
1252 // Apply given scaling factor
1253 if (is_CC) { xsec *= fXSecScaleCC; }
1254 else if (is_NC) { xsec *= fXSecScaleNC; }
1255
1256 // If requested return the free nucleon xsec even for input nuclear tgt
1257 if ( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
1258
1259 int Z = target.Z();
1260 int A = target.A();
1261 int N = A-Z;
1262
1263 // Take into account the number of scattering centers in the target
1264 int NNucl = (SppChannel::InitStateNucleon(spp_channel) == kPdgProton) ? Z : N;
1265 xsec*=NNucl; // nuclear xsec (no nuclear suppression symmetry_factor)
1266
1267 if ( fUsePauliBlocking && A!=1 && kps == kPSWQ2ctpfE )
1268 {
1269 // Calculation of Pauli blocking according to refs. 9-11
1270 double P_Fermi = 0;
1271
1272 // Maximum value of Fermi momentum of target nucleon (GeV)
1273 if ( A<6 || ! fUseRFGParametrization )
1274 {
1275 // look up the Fermi momentum for this target
1276 FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
1277 const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
1278 P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucpdgc);
1279 }
1280 else {
1281 // define the Fermi momentum for this target
1282 P_Fermi = utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(target);
1283 // correct the Fermi momentum for the struck nucleon
1284 if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); }
1285 else { P_Fermi *= TMath::Power( 2.*N/A, 1./3); }
1286 }
1287
1288 double FactorPauli_RES = 1;
1289
1290 if (P_Fermi > 0)
1291 {
1292 if ( 2*P_Fermi < abs_mom_k-abs_mom_q )
1293 FactorPauli_RES = 1;
1294 if ( 2*P_Fermi >= abs_mom_k+abs_mom_q )
1295 FactorPauli_RES = ((3*abs_mom_k*abs_mom_k + abs_mom_q*abs_mom_q)/(2*P_Fermi) - (5*TMath::Power(abs_mom_k,4) + TMath::Power(abs_mom_q,4) + 10*abs_mom_k*abs_mom_k*abs_mom_q*abs_mom_q)/(40*TMath::Power(P_Fermi,3)))/(2*abs_mom_k);
1296 if ( 2*P_Fermi >= abs_mom_k-abs_mom_q && 2*P_Fermi <= abs_mom_k+abs_mom_q )
1297 FactorPauli_RES = ((abs_mom_q + abs_mom_k)*(abs_mom_q + abs_mom_k) - 4*P_Fermi*P_Fermi/5 - TMath::Power(abs_mom_k - abs_mom_q, 3)/(2*P_Fermi)+TMath::Power(abs_mom_k - abs_mom_q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*abs_mom_q*abs_mom_k);
1298 }
1299 xsec *= FactorPauli_RES;
1300 }
1301 return xsec;
1302
1303}
1304
1305//____________________________________________________________________________
1306double MKSPPPXSec2020::Integral(const Interaction * interaction) const
1307{
1308 double xsec = fXSecIntegrator->Integrate(this,interaction);
1309 return xsec;
1310}
1311//____________________________________________________________________________
1312bool MKSPPPXSec2020::ValidProcess(const Interaction * interaction) const
1313{
1314
1315 if(interaction->TestBit(kISkipProcessChk)) return true;
1316
1317 //-- Get the requested SPP channel
1318 SppChannel_t spp_channel = SppChannel::FromInteraction(interaction);
1319 if( spp_channel == kSppNull ) {
1320 return false;
1321 }
1322
1323 return true;
1324
1325}
1326//____________________________________________________________________________
1327inline int MKSPPPXSec2020::Lambda (NucleonPolarization l) const
1328{
1329 return 2*l - 1;
1330}
1331//____________________________________________________________________________
1332inline int MKSPPPXSec2020::Lambda (BosonPolarization l) const
1333{
1334 return (l*(l*(4*l-21)+29)-6)/3;
1335}
1336//____________________________________________________________________________
1337inline int MKSPPPXSec2020::PhaseFactor(BosonPolarization lk, NucleonPolarization l1, NucleonPolarization l2) const
1338{
1339 return lk*(lk*(4*lk - 21) + 29)/6 - l1 - l2;
1340}
1341//____________________________________________________________________________
1342bool MKSPPPXSec2020::ValidKinematics(const Interaction * interaction) const
1343{
1344 // call only after ValidProcess
1345 if ( interaction->TestBit(kISkipKinematicChk) ) return true;
1346
1347 const KPhaseSpace& kps = interaction->PhaseSpace();
1348
1349 // Get kinematical parameters
1350 const InitialState & init_state = interaction -> InitState();
1351 const Kinematics & kinematics = interaction -> Kine();
1352 double Enu = init_state.ProbeE(kRfHitNucRest);
1353 double W = kinematics.W();
1354 double Q2 = kinematics.Q2();
1355
1356 if (Enu < kps.Threshold_SPP_iso())
1357 return false;
1358
1359 Range1D_t Wl = kps.WLim_SPP_iso();
1360 if (fWmax >= Wl.min)
1361 Wl.max = TMath::Min(fWmax, Wl.max);
1362 Range1D_t Q2l = kps.Q2Lim_W_SPP_iso();
1363
1364 if (W < Wl.min || W > Wl.max)
1365 return false;
1366
1367 if (Q2 < Q2l.min || Q2 > Q2l.max)
1368 return false;
1369
1370 return true;
1371
1372}
1373//____________________________________________________________________________
1374void MKSPPPXSec2020::Configure(const Registry & config)
1375{
1376 Algorithm::Configure(config);
1377 this->LoadConfig();
1378}
1379//____________________________________________________________________________
1380void MKSPPPXSec2020::Configure(string config)
1381{
1382 Algorithm::Configure(config);
1383 this->LoadConfig();
1384}
1385//____________________________________________________________________________
1386void MKSPPPXSec2020::LoadConfig(void)
1387{
1388
1389 fResList.Clear();
1390 string resonances ;
1391 GetParam( "ResonanceNameList", resonances ) ;
1392 fResList.DecodeFromNameList(resonances);
1393
1394 // Cross section scaling symmetry_factors
1395 this->GetParam( "RES-CC-XSecScale", fXSecScaleCC );
1396 this->GetParam( "RES-NC-XSecScale", fXSecScaleNC );
1397
1398 // Load all configuration data or set defaults
1399 this->GetParamDef( "RES-Omega" , fOmega, 1.05);
1400
1401 double ma, mv;
1402 this->GetParam( "RES-Ma" , ma);
1403 this->GetParam( "RES-Mv" , mv);
1404 fMa2 = ma*ma;
1405 fMv2 = mv*mv;
1406 this->GetParam( "RES-CA50" , fCA50) ;
1407
1408 this->GetParam( "GVCAL-Cv3" , fCv3) ;
1409 this->GetParam( "GVCAL-Cv4" , fCv4) ;
1410 this->GetParam( "GVCAL-Cv51" , fCv51) ;
1411 this->GetParam( "GVCAL-Cv52" , fCv52) ;
1412
1413 this->GetParam( "FermiConstant", fFermiConstant );
1414
1415 double thw;
1416 this->GetParam( "WeinbergAngle", thw );
1417 fSin2Wein = TMath::Power( TMath::Sin(thw), 2 );
1418
1419 this->GetParam("CKM-Vud", fVud );
1420
1421 // Load all the sub-algorithms needed
1422 fHAmplModelCC = 0;
1423 fHAmplModelNCp = 0;
1424 fHAmplModelNCn = 0;
1425
1426 fHAmplModelCC = dynamic_cast<const RSHelicityAmplModelI *> ( this -> SubAlg( "HelictyAmplCCAlg" ) );
1427 fHAmplModelNCp = dynamic_cast<const RSHelicityAmplModelI *> ( this -> SubAlg( "HelictyAmplNCpAlg" ));
1428 fHAmplModelNCn = dynamic_cast<const RSHelicityAmplModelI *> ( this -> SubAlg( "HelictyAmplNCnAlg" ));
1429
1430 assert(fHAmplModelCC);
1431 assert(fHAmplModelNCp);
1432 assert(fHAmplModelNCn);
1433
1434 // Parameters for calculation of background contribution
1435 fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (this->SubAlg("CCFormFactorsAlg"));
1436 assert(fFormFactorsModel);
1437 fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
1438
1439 fEMFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (this->SubAlg("EMFormFactorsAlg"));
1440 assert(fEMFormFactorsModel);
1441 fEMFormFactors.SetModel(fEMFormFactorsModel); // <-- attach algorithm
1442
1443 this->GetParam("FermiMomentumTable", fKFTable);
1444 this->GetParam("RFG-UseParametrization", fUseRFGParametrization);
1445 this->GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
1446
1447
1448 this->GetParamDef("MK-fPi", f_pi, 0.093 );
1449 this->GetParamDef("QEL-FA0", FA0, -1.26 );
1450 this->GetParamDef("MK-Frho0", Frho0, 1.0 );
1451
1452 this->GetParamDef("MK-NonResBkg-VWmin", fBkgVWmin, 1.30 );
1453 this->GetParamDef("MK-NonResBkg-VWmax", fBkgVWmax, 1.60 );
1454 this->GetParamDef("MK-NonResBkg-V3", fBkgV3, 8.08497 );
1455 this->GetParamDef("MK-NonResBkg-V2", fBkgV2, -41.6767 );
1456 this->GetParamDef("MK-NonResBkg-V1", fBkgV1, 66.3419 );
1457 this->GetParamDef("MK-NonResBkg-V0", fBkgV0, -32.5733 );
1458 this->GetParamDef("Mass-rho770", fRho770Mass, 0.7758 );
1459 this->GetParamDef("MK-WMax", fWmax, -1.0 );
1460
1461 this->GetParamDef("UseAuthorCode", fUseAuthorCode, false );
1462
1463 // Load the differential cross section integrator
1464 fXSecIntegrator = dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
1465 assert(fXSecIntegrator);
1466
1467}
1468//____________________________________________________________________________
pFATAL
#define pFATAL
Definition Messenger.h:56
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils.
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition Algorithm.cxx:345
genie::FermiMomentumTablePool
Singleton class to load & serve tables of Fermi momentum constants.
Definition FermiMomentumTablePool.h:35
genie::FermiMomentumTablePool::GetTable
const FermiMomentumTable * GetTable(string name)
Definition FermiMomentumTablePool.cxx:76
genie::FermiMomentumTablePool::Instance
static FermiMomentumTablePool * Instance(void)
Definition FermiMomentumTablePool.cxx:62
genie::FermiMomentumTable
A table of Fermi momentum constants.
Definition FermiMomentumTable.h:34
genie::FermiMomentumTable::FindClosestKF
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
Definition FermiMomentumTable.cxx:45
genie::InitialState
Initial State information.
Definition InitialState.h:48
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition InitialState.h:66
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition InitialState.h:64
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition InitialState.cxx:384
genie::Interaction
Summary information for an interaction.
Definition Interaction.h:56
genie::Interaction::FSPrimLepton
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Definition Interaction.cxx:126
genie::Interaction::PhaseSpace
const KPhaseSpace & PhaseSpace(void) const
Definition Interaction.h:73
genie::KPhaseSpace
Kinematical phase space.
Definition KPhaseSpace.h:33
genie::KPhaseSpace::Threshold_SPP_iso
double Threshold_SPP_iso(void) const
Energy limit for resonance single pion production on isoscalar nucleon.
Definition KPhaseSpace.cxx:1000
genie::KPhaseSpace::WLim_SPP_iso
Range1D_t WLim_SPP_iso(void) const
W limits for resonance single pion production on isoscalar nucleon.
Definition KPhaseSpace.cxx:1043
genie::KPhaseSpace::Q2Lim_W_SPP_iso
Range1D_t Q2Lim_W_SPP_iso(void) const
Q2 limits @ fixed W for resonance single pion production on isoscalar nucleon.
Definition KPhaseSpace.cxx:1111
genie::KinePhaseSpace::AsString
static string AsString(KinePhaseSpace_t kps)
Definition KinePhaseSpace.h:85
genie::Kinematics
Generated/set kinematical variables for an event.
Definition Kinematics.h:39
genie::Kinematics::Q2
double Q2(bool selected=false) const
Definition Kinematics.cxx:125
genie::Kinematics::GetKV
double GetKV(KineVar_t kv) const
Definition Kinematics.cxx:323
genie::Kinematics::W
double W(bool selected=false) const
Definition Kinematics.cxx:157
genie::MKSPPPXSec2020::fRho770Mass
double fRho770Mass
Mass of rho(770) meson.
Definition MKSPPPXSec2020.h:412
genie::MKSPPPXSec2020::Integral
double Integral(const Interaction *i) const
Definition MKSPPPXSec2020.cxx:1306
genie::MKSPPPXSec2020::fHAmplModelNCp
const RSHelicityAmplModelI * fHAmplModelNCp
Definition MKSPPPXSec2020.h:370
genie::MKSPPPXSec2020::Configure
void Configure(const Registry &config)
Definition MKSPPPXSec2020.cxx:1374
genie::MKSPPPXSec2020::fSin2Wein
double fSin2Wein
sin^2(Weingberg angle)
Definition MKSPPPXSec2020.h:384
genie::MKSPPPXSec2020::fUseRFGParametrization
bool fUseRFGParametrization
Use parametrization for fermi momentum insted of table?
Definition MKSPPPXSec2020.h:393
genie::MKSPPPXSec2020::fHAmplModelNCn
const RSHelicityAmplModelI * fHAmplModelNCn
Definition MKSPPPXSec2020.h:371
genie::MKSPPPXSec2020::fUsePauliBlocking
bool fUsePauliBlocking
Account for Pauli blocking?
Definition MKSPPPXSec2020.h:394
genie::MKSPPPXSec2020::fCA50
double fCA50
FKR parameter Zeta.
Definition MKSPPPXSec2020.h:376
genie::MKSPPPXSec2020::fMa2
double fMa2
(axial mass)^2
Definition MKSPPPXSec2020.h:378
genie::MKSPPPXSec2020::fMv2
double fMv2
(vector mass)^2
Definition MKSPPPXSec2020.h:379
genie::MKSPPPXSec2020::fHAmplModelCC
const RSHelicityAmplModelI * fHAmplModelCC
Definition MKSPPPXSec2020.h:369
genie::MKSPPPXSec2020::fFormFactorsModel
const QELFormFactorsModelI * fFormFactorsModel
Quasielastic form facors model for background contribution.
Definition MKSPPPXSec2020.h:389
genie::MKSPPPXSec2020::ValidKinematics
bool ValidKinematics(const Interaction *interaction) const
Is the input kinematical point a physically allowed one?
Definition MKSPPPXSec2020.cxx:1342
genie::MKSPPPXSec2020::fEMFormFactorsModel
const QELFormFactorsModelI * fEMFormFactorsModel
Electromagnetic form factors model for background contribution.
Definition MKSPPPXSec2020.h:390
genie::MKSPPPXSec2020::fWmax
double fWmax
The value above which the partial cross section is set to zero (if negative then not appliciable)
Definition MKSPPPXSec2020.h:413
genie::MKSPPPXSec2020::fXSecScaleNC
double fXSecScaleNC
External NC xsec scaling factor.
Definition MKSPPPXSec2020.h:387
genie::MKSPPPXSec2020::NucleonPolarizationIterator
Iterator< NucleonPolarization, NucleonPolarization::MINUS, NucleonPolarization::PLUS > NucleonPolarizationIterator
Definition MKSPPPXSec2020.h:155
genie::MKSPPPXSec2020::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition MKSPPPXSec2020.h:417
genie::MKSPPPXSec2020::fEMFormFactors
QELFormFactors fEMFormFactors
Electromagnetic form facors for background contribution.
Definition MKSPPPXSec2020.h:397
genie::MKSPPPXSec2020::fVud
double fVud
|Vud| (magnitude ud-element of CKM-matrix)
Definition MKSPPPXSec2020.h:385
genie::MKSPPPXSec2020::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition MKSPPPXSec2020.cxx:1312
genie::MKSPPPXSec2020::fCv3
double fCv3
GV calculation coeffs.
Definition MKSPPPXSec2020.h:380
genie::MKSPPPXSec2020::fXSecScaleCC
double fXSecScaleCC
External CC xsec scaling factor.
Definition MKSPPPXSec2020.h:386
genie::MKSPPPXSec2020::FA0
double FA0
Axial coupling (value of axial form factor at Q2=0)
Definition MKSPPPXSec2020.h:399
genie::MKSPPPXSec2020::fKFTable
string fKFTable
Table of Fermi momentum (kF) constants for various nuclei.
Definition MKSPPPXSec2020.h:392
genie::MKSPPPXSec2020::fUseAuthorCode
bool fUseAuthorCode
Use author code?
Definition MKSPPPXSec2020.h:415
genie::MKSPPPXSec2020::f_pi
double f_pi
Constant for pion-nucleon interaction.
Definition MKSPPPXSec2020.h:398
genie::MKSPPPXSec2020::fOmega
double fOmega
FKR parameter Omega.
Definition MKSPPPXSec2020.h:377
genie::MKSPPPXSec2020::fFormFactors
QELFormFactors fFormFactors
Quasielastic form facors for background contribution.
Definition MKSPPPXSec2020.h:396
genie::MKSPPPXSec2020::Lambda
int Lambda(NucleonPolarization l) const
Definition MKSPPPXSec2020.cxx:1327
genie::MKSPPPXSec2020::BosonPolarizationIterator
Iterator< BosonPolarization, BosonPolarization::LEFT, BosonPolarization::PLUS0 > BosonPolarizationIterator
Definition MKSPPPXSec2020.h:154
genie::MKSPPPXSec2020::PhaseFactor
int PhaseFactor(BosonPolarization lk, NucleonPolarization l1, NucleonPolarization l2) const
Definition MKSPPPXSec2020.cxx:1337
genie::MKSPPPXSec2020::XSec
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition MKSPPPXSec2020.cxx:69
genie::PDGLibrary
Singleton class to load & serve a TDatabasePDG.
Definition PDGLibrary.h:36
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition PDGLibrary.cxx:68
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition PDGLibrary.cxx:86
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition ProcessInfo.h:46
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition ProcessInfo.cxx:208
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition ProcessInfo.cxx:203
genie::QELFormFactorsModelI
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
Definition QELFormFactorsModelI.h:28
genie::RSHelicityAmplModelI
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
Definition RSHelicityAmplModelI.h:28
genie::RSHelicityAmplModelI::Compute
virtual const RSHelicityAmpl & Compute(Resonance_t res, const FKR &fkr) const =0
genie::RSHelicityAmpl
A class holding the Rein-Sehgal's helicity amplitudes.
Definition RSHelicityAmpl.h:40
genie::RSHelicityAmpl::Amp0Minus
double Amp0Minus(void) const
Definition RSHelicityAmpl.h:59
genie::RSHelicityAmpl::AmpPlus3
double AmpPlus3(void) const
Definition RSHelicityAmpl.h:58
genie::RSHelicityAmpl::Amp0Plus
double Amp0Plus(void) const
Definition RSHelicityAmpl.h:60
genie::RSHelicityAmpl::AmpPlus1
double AmpPlus1(void) const
Definition RSHelicityAmpl.h:56
genie::RSHelicityAmpl::AmpMinus3
double AmpMinus3(void) const
Definition RSHelicityAmpl.h:57
genie::RSHelicityAmpl::AmpMinus1
double AmpMinus1(void) const
return helicity amplitude
Definition RSHelicityAmpl.h:55
genie::Range1D_t
A simple [min,max] interval for doubles.
Definition Range1.h:43
genie::Range1D_t::max
double max
Definition Range1.h:53
genie::Range1D_t::min
double min
Definition Range1.h:52
genie::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition Registry.h:65
genie::SppChannel::Isospin1Coefficients
static double Isospin1Coefficients(SppChannel_t channel)
Definition SppChannel.h:317
genie::SppChannel::Isospin3Coefficients
static double Isospin3Coefficients(SppChannel_t channel)
Definition SppChannel.h:280
genie::SppChannel::FromInteraction
static SppChannel_t FromInteraction(const Interaction *interaction)
Definition SppChannel.h:402
genie::SppChannel::FinStateIsospin
static int FinStateIsospin(SppChannel_t channel)
Definition SppChannel.h:211
genie::SppChannel::InitStateNucleon
static int InitStateNucleon(SppChannel_t channel)
Definition SppChannel.h:103
genie::SppChannel::BranchingRatio
static double BranchingRatio(Resonance_t res)
Definition SppChannel.h:357
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition Target.h:40
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition Target.cxx:304
genie::Target::Z
int Z(void) const
Definition Target.h:68
genie::Target::A
int A(void) const
Definition Target.h:70
genie::XSecIntegratorI
Cross Section Integrator Interface.
Definition XSecIntegratorI.h:29
genie::constants
Basic constants.
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition PDGUtils.cxx:110
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition PDGUtils.cxx:71
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition PDGUtils.cxx:336
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition PDGUtils.cxx:118
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition KineUtils.cxx:130
genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
Definition NuclearUtils.cxx:501
genie::utils::res::AngularMom
int AngularMom(Resonance_t res)
Definition BaryonResUtils.cxx:622
genie::utils::res::Cjsgn_plus
int Cjsgn_plus(Resonance_t res)
Definition BaryonResUtils.cxx:652
genie::utils::res::Width
double Width(Resonance_t res)
resonance width (GeV)
Definition BaryonResUtils.cxx:463
genie::utils::res::Dsgn
int Dsgn(Resonance_t res)
Definition BaryonResUtils.cxx:685
genie::utils::res::Isospin
int Isospin(Resonance_t res)
Definition BaryonResUtils.cxx:591
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition BaryonResUtils.cxx:436
genie::utils::res::OrbitalAngularMom
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
Definition BaryonResUtils.cxx:531
genie::utils::res::ResonanceIndex
int ResonanceIndex(Resonance_t res)
resonance idx, quark model / SU(6)
Definition BaryonResUtils.cxx:561
genie
THE MAIN GENIE PROJECT NAMESPACE
Definition AlgCmp.h:25
genie::kPdgPiM
const int kPdgPiM
Definition PDGCodes.h:159
genie::gAbortingInErr
bool gAbortingInErr
Definition Messenger.cxx:34
genie::kPdgProton
const int kPdgProton
Definition PDGCodes.h:81
genie::kPdgPi0
const int kPdgPi0
Definition PDGCodes.h:160
genie::kPSWQ2ctpphipfE
@ kPSWQ2ctpphipfE
Definition KinePhaseSpace.h:76
genie::kPdgNeutron
const int kPdgNeutron
Definition PDGCodes.h:83
genie::kKVphip
@ kKVphip
Definition KineVar.h:54
genie::kKVctp
@ kKVctp
Definition KineVar.h:55
genie::kSpp_vp_cc_10100
@ kSpp_vp_cc_10100
Definition SppChannel.h:50
genie::kSpp_vbn_nc_10001
@ kSpp_vbn_nc_10001
Definition SppChannel.h:66
genie::kSppNull
@ kSppNull
Definition SppChannel.h:46
genie::kSpp_vbp_nc_10010
@ kSpp_vbp_nc_10010
Definition SppChannel.h:63
genie::kSpp_vp_nc_01100
@ kSpp_vp_nc_01100
Definition SppChannel.h:55
genie::kSpp_vn_nc_01010
@ kSpp_vn_nc_01010
Definition SppChannel.h:56
genie::kSpp_vbn_cc_01001
@ kSpp_vbn_cc_01001
Definition SppChannel.h:59
genie::kSpp_vn_cc_01100
@ kSpp_vn_cc_01100
Definition SppChannel.h:52
genie::kSpp_vn_nc_10001
@ kSpp_vn_nc_10001
Definition SppChannel.h:57
genie::kSpp_vbp_nc_01100
@ kSpp_vbp_nc_01100
Definition SppChannel.h:64
genie::kSpp_vn_cc_10010
@ kSpp_vn_cc_10010
Definition SppChannel.h:51
genie::kSpp_vbp_cc_01010
@ kSpp_vbp_cc_01010
Definition SppChannel.h:60
genie::kSpp_vbp_cc_10001
@ kSpp_vbp_cc_10001
Definition SppChannel.h:61
genie::kSpp_vbn_nc_01010
@ kSpp_vbn_nc_01010
Definition SppChannel.h:65
genie::kSpp_vp_nc_10010
@ kSpp_vp_nc_10010
Definition SppChannel.h:54
genie::KinePhaseSpace_t
enum genie::EKinePhaseSpace KinePhaseSpace_t
genie::kISkipKinematicChk
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition Interaction.h:48
genie::SppChannel_t
enum genie::ESppChannel SppChannel_t
genie::kRfHitNucRest
@ kRfHitNucRest
Definition RefFrame.h:30
genie::kPdgPiP
const int kPdgPiP
Definition PDGCodes.h:158
genie::kISkipProcessChk
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition Interaction.h:47
genie::kIAssumeFreeNucleon
const UInt_t kIAssumeFreeNucleon
Definition Interaction.h:49
Generated on for GENIEGenerator by doxygen 1.14.0