GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
MKSPPPXSec2020.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, 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 
48 using namespace genie;
49 using namespace genie::constants;
50 
51 //____________________________________________________________________________
52 MKSPPPXSec2020::MKSPPPXSec2020() :
53 XSecAlgorithmI("genie::MKSPPPXSec2020")
54 {
55 
56 }
57 //____________________________________________________________________________
58 MKSPPPXSec2020::MKSPPPXSec2020(string config) :
59 XSecAlgorithmI("genie::MKSPPPXSec2020", config)
60 {
61 
62 }
63 //____________________________________________________________________________
64 MKSPPPXSec2020::~MKSPPPXSec2020()
65 {
66 
67 }
68 //____________________________________________________________________________
69 double 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 //____________________________________________________________________________
1306 double MKSPPPXSec2020::Integral(const Interaction * interaction) const
1307 {
1308  double xsec = fXSecIntegrator->Integrate(this,interaction);
1309  return xsec;
1310 }
1311 //____________________________________________________________________________
1312 bool 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 //____________________________________________________________________________
1327 inline int MKSPPPXSec2020::Lambda (NucleonPolarization l) const
1328 {
1329  return 2*l - 1;
1330 }
1331 //____________________________________________________________________________
1332 inline int MKSPPPXSec2020::Lambda (BosonPolarization l) const
1333 {
1334  return (l*(l*(4*l-21)+29)-6)/3;
1335 }
1336 //____________________________________________________________________________
1337 inline int MKSPPPXSec2020::PhaseFactor(BosonPolarization lk, NucleonPolarization l1, NucleonPolarization l2) const
1338 {
1339  return lk*(lk*(4*lk - 21) + 29)/6 - l1 - l2;
1340 }
1341 //____________________________________________________________________________
1342 bool 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 //____________________________________________________________________________
1374 void MKSPPPXSec2020::Configure(const Registry & config)
1375 {
1376  Algorithm::Configure(config);
1377  this->LoadConfig();
1378 }
1379 //____________________________________________________________________________
1380 void MKSPPPXSec2020::Configure(string config)
1381 {
1382  Algorithm::Configure(config);
1383  this->LoadConfig();
1384 }
1385 //____________________________________________________________________________
1386 void 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 //____________________________________________________________________________
genie::XSecAlgorithmI
Cross Section Calculation Interface.
Definition: XSecAlgorithmI.h:27
genie::SppChannel::FromInteraction
static SppChannel_t FromInteraction(const Interaction *interaction)
Definition: SppChannel.h:402
genie::Kinematics::W
double W(bool selected=false) const
Definition: Kinematics.cxx:157
genie::RSHelicityAmpl::AmpMinus3
double AmpMinus3(void) const
Definition: RSHelicityAmpl.h:57
genie::Interaction::PhaseSpace
const KPhaseSpace & PhaseSpace(void) const
Definition: Interaction.h:73
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:203
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
genie::utils::mec::J
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
genie::units::g
static constexpr double g
Definition: Units.h:144
genie::FKR::Rminus
double Rminus
Definition: FKR.h:50
genie::MKSPPPXSec2020::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition: MKSPPPXSec2020.h:417
genie::XSecIntegratorI
Cross Section Integrator Interface.
Definition: XSecIntegratorI.h:29
genie::utils::kinematics::Q2
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
genie::Target::HitNucPdg
int HitNucPdg(void) const
Definition: Target.cxx:304
genie::MKSPPPXSec2020::BosonPolarizationIterator
Iterator< BosonPolarization, BosonPolarization::LEFT, BosonPolarization::PLUS0 > BosonPolarizationIterator
Definition: MKSPPPXSec2020.h:154
genie::FKR::Ra
double Ra
Definition: FKR.h:42
genie::Target::A
int A(void) const
Definition: Target.h:70
genie::Range1D_t
A simple [min,max] interval for doubles.
Definition: Range1.h:42
genie::QELFormFactors::SetModel
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
Definition: QELFormFactors.cxx:43
genie::RSHelicityAmplModelI::Compute
virtual const RSHelicityAmpl & Compute(Resonance_t res, const FKR &fkr) const =0
genie::RSHelicityAmpl::AmpPlus3
double AmpPlus3(void) const
Definition: RSHelicityAmpl.h:58
pFATAL
#define pFATAL
Definition: Messenger.h:56
genie::BaryonResList::DecodeFromNameList
void DecodeFromNameList(string list, string delimiter=",")
Definition: BaryonResList.cxx:101
genie::units::s
static constexpr double s
Definition: Units.h:95
genie::SppChannel::Isospin1Coefficients
static double Isospin1Coefficients(SppChannel_t channel)
Definition: SppChannel.h:317
genie::FermiMomentumTablePool::Instance
static FermiMomentumTablePool * Instance(void)
Definition: FermiMomentumTablePool.cxx:62
genie::Kinematics
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
genie::MKSPPPXSec2020::fSin2Wein
double fSin2Wein
sin^2(Weingberg angle)
Definition: MKSPPPXSec2020.h:384
genie::FKR::Lamda
double Lamda
Definition: FKR.h:37
genie::utils::res::Mass
double Mass(Resonance_t res)
resonance mass (GeV)
Definition: BaryonResUtils.cxx:436
genie::FKR::R
double R
Definition: FKR.h:45
genie::FermiMomentumTable
A table of Fermi momentum constants.
Definition: FermiMomentumTable.h:33
genie::utils::res::Width
double Width(Resonance_t res)
resonance width (GeV)
Definition: BaryonResUtils.cxx:463
genie::MKSPPPXSec2020::fHAmplModelCC
const RSHelicityAmplModelI * fHAmplModelCC
Definition: MKSPPPXSec2020.h:369
genie::KPhaseSpace::Threshold_SPP_iso
double Threshold_SPP_iso(void) const
Energy limit for resonance single pion production on isoscalar nucleon.
Definition: KPhaseSpace.cxx:999
genie::utils::res::AngularMom
int AngularMom(Resonance_t res)
Definition: BaryonResUtils.cxx:622
genie::KinePhaseSpace_t
enum genie::EKinePhaseSpace KinePhaseSpace_t
genie::MKSPPPXSec2020::fUseRFGParametrization
bool fUseRFGParametrization
Use parametrization for fermi momentum insted of table?
Definition: MKSPPPXSec2020.h:393
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::fMa2
double fMa2
(axial mass)^2
Definition: MKSPPPXSec2020.h:378
genie::MKSPPPXSec2020::fXSecScaleNC
double fXSecScaleNC
External NC xsec scaling factor.
Definition: MKSPPPXSec2020.h:387
genie::MKSPPPXSec2020::Integral
double Integral(const Interaction *i) const
Definition: MKSPPPXSec2020.cxx:1306
genie::MKSPPPXSec2020::fCA50
double fCA50
FKR parameter Zeta.
Definition: MKSPPPXSec2020.h:376
genie::SppChannel_t
enum genie::ESppChannel SppChannel_t
genie::utils::res::Cjsgn_plus
int Cjsgn_plus(Resonance_t res)
Definition: BaryonResUtils.cxx:652
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
genie::RSHelicityAmpl::Amp0Plus
double Amp0Plus(void) const
Definition: RSHelicityAmpl.h:60
genie::units::fm3
static constexpr double fm3
Definition: Units.h:77
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
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:1110
genie::FKR::Tv
double Tv
Definition: FKR.h:38
genie::MKSPPPXSec2020::fResList
BaryonResList fResList
Definition: MKSPPPXSec2020.h:419
genie::RSHelicityAmpl
A class holding the Rein-Sehgal&#39;s helicity amplitudes.
Definition: RSHelicityAmpl.h:40
genie::pdg::IsProton
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
genie::ProcessInfo::IsWeakNC
bool IsWeakNC(void) const
Definition: ProcessInfo.cxx:208
genie::FermiMomentumTablePool
Singleton class to load &amp; serve tables of Fermi momentum constants.
Definition: FermiMomentumTablePool.h:34
LOG
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
genie::QELFormFactorsModelI
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
Definition: QELFormFactorsModelI.h:28
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::SppChannel::InitStateNucleon
static int InitStateNucleon(SppChannel_t channel)
Definition: SppChannel.h:103
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::FermiMomentumTablePool::GetTable
const FermiMomentumTable * GetTable(string name)
Definition: FermiMomentumTablePool.cxx:76
genie::KinePhaseSpace::AsString
static string AsString(KinePhaseSpace_t kps)
Definition: KinePhaseSpace.h:85
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
genie::MKSPPPXSec2020::Lambda
int Lambda(NucleonPolarization l) const
Definition: MKSPPPXSec2020.cxx:1327
genie::FKR::T
double T
Definition: FKR.h:46
genie::FKR::Rv
double Rv
Definition: FKR.h:39
genie::MKSPPPXSec2020::fHAmplModelNCp
const RSHelicityAmplModelI * fHAmplModelNCp
Definition: MKSPPPXSec2020.h:370
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
genie::Target
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
genie::Algorithm::Configure
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
genie::InitialState::ProbePdg
int ProbePdg(void) const
Definition: InitialState.h:64
genie::KPhaseSpace
Kinematical phase space.
Definition: KPhaseSpace.h:33
genie::MKSPPPXSec2020::fRho770Mass
double fRho770Mass
Mass of rho(770) meson.
Definition: MKSPPPXSec2020.h:412
genie::RSHelicityAmpl::AmpMinus1
double AmpMinus1(void) const
return helicity amplitude
Definition: RSHelicityAmpl.h:55
genie::kPdgPiP
const int kPdgPiP
Definition: PDGCodes.h:158
genie::kPdgPi0
const int kPdgPi0
Definition: PDGCodes.h:160
genie::MKSPPPXSec2020::fMv2
double fMv2
(vector mass)^2
Definition: MKSPPPXSec2020.h:379
genie::utils::res::OrbitalAngularMom
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
Definition: BaryonResUtils.cxx:531
genie::RSHelicityAmpl::Amp0Minus
double Amp0Minus(void) const
Definition: RSHelicityAmpl.h:59
genie::MKSPPPXSec2020::fHAmplModelNCn
const RSHelicityAmplModelI * fHAmplModelNCn
Definition: MKSPPPXSec2020.h:371
genie::Target::Z
int Z(void) const
Definition: Target.h:68
genie::MKSPPPXSec2020::FA0
double FA0
Axial coupling (value of axial form factor at Q2=0)
Definition: MKSPPPXSec2020.h:399
genie::Kinematics::GetKV
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
genie::MKSPPPXSec2020::fEMFormFactorsModel
const QELFormFactorsModelI * fEMFormFactorsModel
Electromagnetic form factors model for background contribution.
Definition: MKSPPPXSec2020.h:390
genie::MKSPPPXSec2020::f_pi
double f_pi
Constant for pion-nucleon interaction.
Definition: MKSPPPXSec2020.h:398
genie::RSHelicityAmplModelI
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
Definition: RSHelicityAmplModelI.h:27
genie::QELFormFactors::xiF2V
double xiF2V(void) const
Get the computed form factor xi*F2V.
Definition: QELFormFactors.h:56
genie::kISkipKinematicChk
const UInt_t kISkipKinematicChk
if set, skip kinematic validity checks
Definition: Interaction.h:48
genie::MKSPPPXSec2020::Configure
void Configure(const Registry &config)
Definition: MKSPPPXSec2020.cxx:1374
genie::Interaction::FSPrimLepton
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Definition: Interaction.cxx:126
genie::QELFormFactors::Calculate
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
Definition: QELFormFactors.cxx:49
genie::utils::res::Dsgn
int Dsgn(Resonance_t res)
Definition: BaryonResUtils.cxx:685
genie::SppChannel::Isospin3Coefficients
static double Isospin3Coefficients(SppChannel_t channel)
Definition: SppChannel.h:280
genie::MKSPPPXSec2020::fKFTable
string fKFTable
Table of Fermi momentum (kF) constants for various nuclei.
Definition: MKSPPPXSec2020.h:392
genie::FKR::C
double C
Definition: FKR.h:44
genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
Definition: NuclearUtils.cxx:501
genie::Range1D_t::max
double max
Definition: Range1.h:53
genie::FKR::Tplus
double Tplus
Definition: FKR.h:47
genie::PDGLibrary::Instance
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
genie::MKSPPPXSec2020::fCv3
double fCv3
GV calculation coeffs.
Definition: MKSPPPXSec2020.h:380
genie::MKSPPPXSec2020::fUsePauliBlocking
bool fUsePauliBlocking
Account for Pauli blocking?
Definition: MKSPPPXSec2020.h:394
genie::FKR::B
double B
Definition: FKR.h:43
genie::PDGLibrary
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
genie::FKR::Rplus
double Rplus
Definition: FKR.h:49
genie::Registry
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
genie::kIAssumeFreeNucleon
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
genie::MKSPPPXSec2020::fOmega
double fOmega
FKR parameter Omega.
Definition: MKSPPPXSec2020.h:377
genie::MKSPPPXSec2020::fVud
double fVud
|Vud| (magnitude ud-element of CKM-matrix)
Definition: MKSPPPXSec2020.h:385
genie::FKR::Tminus
double Tminus
Definition: FKR.h:48
genie::RSHelicityAmpl::AmpPlus1
double AmpPlus1(void) const
Definition: RSHelicityAmpl.h:56
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
genie::MKSPPPXSec2020::PhaseFactor
int PhaseFactor(BosonPolarization lk, NucleonPolarization l1, NucleonPolarization l2) const
Definition: MKSPPPXSec2020.cxx:1337
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:1042
genie::kPdgPiM
const int kPdgPiM
Definition: PDGCodes.h:159
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
genie::MKSPPPXSec2020::fUseAuthorCode
bool fUseAuthorCode
Use author code?
Definition: MKSPPPXSec2020.h:415
genie::Range1D_t::min
double min
Definition: Range1.h:52
genie::kPdgProton
const int kPdgProton
Definition: PDGCodes.h:81
genie::QELFormFactors::F1V
double F1V(void) const
Get the computed form factor F1V.
Definition: QELFormFactors.h:53
genie::PDGLibrary::Find
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
genie::FermiMomentumTable::FindClosestKF
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
Definition: FermiMomentumTable.cxx:45
genie::Kinematics::Q2
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
genie::MKSPPPXSec2020::fFormFactors
QELFormFactors fFormFactors
Quasielastic form facors for background contribution.
Definition: MKSPPPXSec2020.h:396
genie::Algorithm::GetParamDef
bool GetParamDef(const RgKey &name, T &p, const T &def) const
genie::Algorithm::GetParam
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
genie::InitialState::Tgt
const Target & Tgt(void) const
Definition: InitialState.h:66
genie::MKSPPPXSec2020::fXSecScaleCC
double fXSecScaleCC
External CC xsec scaling factor.
Definition: MKSPPPXSec2020.h:386
genie::MKSPPPXSec2020::fEMFormFactors
QELFormFactors fEMFormFactors
Electromagnetic form facors for background contribution.
Definition: MKSPPPXSec2020.h:397
genie::MKSPPPXSec2020::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: MKSPPPXSec2020.cxx:1312
genie::SppChannel::BranchingRatio
static double BranchingRatio(Resonance_t res)
Definition: SppChannel.h:357
genie::XSecIntegratorI::Integrate
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
genie::SppChannel::FinStateIsospin
static int FinStateIsospin(SppChannel_t channel)
Definition: SppChannel.h:211
genie::gAbortingInErr
bool gAbortingInErr
Definition: Messenger.cxx:34
genie::MKSPPPXSec2020::fFormFactorsModel
const QELFormFactorsModelI * fFormFactorsModel
Quasielastic form facors model for background contribution.
Definition: MKSPPPXSec2020.h:389
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition: InitialState.cxx:384
genie::kPdgNeutron
const int kPdgNeutron
Definition: PDGCodes.h:83
genie::QELFormFactors::FA
double FA(void) const
Get the computed form factor FA.
Definition: QELFormFactors.h:59
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
genie::MKSPPPXSec2020::ValidKinematics
bool ValidKinematics(const Interaction *interaction) const
Is the input kinematical point a physically allowed one?
Definition: MKSPPPXSec2020.cxx:1342
genie::FKR::S
double S
Definition: FKR.h:40
genie::FKR::Ta
double Ta
Definition: FKR.h:41
genie::utils::res::Isospin
int Isospin(Resonance_t res)
Definition: BaryonResUtils.cxx:591
genie::kISkipProcessChk
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
genie::utils::res::ResonanceIndex
int ResonanceIndex(Resonance_t res)
resonance idx, quark model / SU(6)
Definition: BaryonResUtils.cxx:561
genie::InitialState
Initial State information.
Definition: InitialState.h:48
genie::MKSPPPXSec2020::NucleonPolarizationIterator
Iterator< NucleonPolarization, NucleonPolarization::MINUS, NucleonPolarization::PLUS > NucleonPolarizationIterator
Definition: MKSPPPXSec2020.h:155
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
Generated on Wed Jul 17 2024 14:46:43 for GENIEGenerator by   doxygen 1.8.5