GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
BSKLNBaseRESPXSec2014.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 
6  Steve Dytman
7  University of Pittsburgh
8 
9  Jarek Nowak
10  University of Lancaster
11 
12  Gabe Perdue
13  Fermilab
14 
15  Costas Andreopoulos <c.andreopoulos \at cern.ch>
16  University of Liverpool
17 
18  Afroditi Papadopoulou <apapadop \at mit.edu>
19  Massachusetts Institute of Technology
20 
21  Adi Ashkenazi <adishka \at gmail.com>
22  Massachusetts Institute of Technology
23 
24  Igor Kakorin <kakorin@jinr.ru>
25  Joint Institute for Nuclear Research
26 */
27 //____________________________________________________________________________
28 
29 #include <TMath.h>
30 #include <TSystem.h>
31 
32 #include "Framework/Algorithm/AlgFactory.h"
33 #include "Framework/Algorithm/AlgConfigPool.h"
34 #include "Framework/ParticleData/BaryonResUtils.h"
35 #include "Framework/Conventions/GBuild.h"
36 #include "Framework/Conventions/Constants.h"
37 #include "Framework/Conventions/RefFrame.h"
38 #include "Framework/Conventions/KineVar.h"
39 #include "Framework/Conventions/Units.h"
40 #include "Framework/Messenger/Messenger.h"
41 #include "Framework/ParticleData/PDGCodes.h"
42 #include "Framework/ParticleData/PDGUtils.h"
43 #include "Framework/Utils/KineUtils.h"
44 #include "Framework/Numerical/MathUtils.h"
45 #include "Framework/Utils/Range1.h"
46 #include "Framework/Utils/BWFunc.h"
47 #include "Physics/XSectionIntegration/XSecIntegratorI.h"
48 #include "Physics/Resonance/XSection/BSKLNBaseRESPXSec2014.h"
49 #include "Physics/Resonance/XSection/RSHelicityAmplModelI.h"
50 #include "Physics/Resonance/XSection/RSHelicityAmpl.h"
51 #include "Physics/NuclearState/FermiMomentumTablePool.h"
52 #include "Physics/NuclearState/FermiMomentumTable.h"
53 #include "Physics/NuclearState/NuclearUtils.h"
54 
55 using namespace genie;
56 using namespace genie::constants;
57 
58 //____________________________________________________________________________
59 BSKLNBaseRESPXSec2014::BSKLNBaseRESPXSec2014(string name) :
60 XSecAlgorithmI(name)
61 {
62 
63 }
64 //____________________________________________________________________________
65 BSKLNBaseRESPXSec2014::BSKLNBaseRESPXSec2014(string name, string config) :
66 XSecAlgorithmI(name, config)
67 {
68 
69 }
70 //____________________________________________________________________________
71 BSKLNBaseRESPXSec2014::~BSKLNBaseRESPXSec2014()
72 {
73 
74 }
75 //____________________________________________________________________________
76 double BSKLNBaseRESPXSec2014::XSec(
77  const Interaction * interaction, KinePhaseSpace_t kps) const
78 {
79  if(! this -> ValidProcess (interaction) ) return 0.;
80  if(! this -> ValidKinematics (interaction) ) return 0.;
81 
82  const InitialState & init_state = interaction -> InitState();
83  const ProcessInfo & proc_info = interaction -> ProcInfo();
84  const Target & target = init_state.Tgt();
85 
86  // Get kinematical parameters
87  const Kinematics & kinematics = interaction -> Kine();
88  double W = kinematics.W();
89  double q2 = kinematics.q2();
90  double costh = kinematics.FSLeptonP4().CosTheta();
91 
92  // Under the DIS/RES joining scheme, xsec(RES)=0 for W>=Wcut
93  if(fUsingDisResJoin) {
94  if(W>=fWcut) {
95 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
96  LOG("BSKLNBaseRESPXSec2014", pDEBUG)
97  << "RES/DIS Join Scheme: XSec[RES, W=" << W
98  << " >= Wcut=" << fWcut << "] = 0";
99 #endif
100  return 0;
101  }
102  }
103 
104  // Get the input baryon resonance
105  Resonance_t resonance = interaction->ExclTag().Resonance();
106  string resname = utils::res::AsString(resonance);
107  bool is_delta = utils::res::IsDelta (resonance);
108 
109  // Get the neutrino, hit nucleon & weak current
110  int nucpdgc = target.HitNucPdg();
111  int probepdgc = init_state.ProbePdg();
112  bool is_nu = pdg::IsNeutrino (probepdgc);
113  bool is_nubar = pdg::IsAntiNeutrino (probepdgc);
114  bool is_lplus = pdg::IsPosChargedLepton (probepdgc);
115  bool is_lminus = pdg::IsNegChargedLepton (probepdgc);
116  bool is_p = pdg::IsProton (nucpdgc);
117  bool is_n = pdg::IsNeutron (nucpdgc);
118  bool is_CC = proc_info.IsWeakCC();
119  bool is_NC = proc_info.IsWeakNC();
120  bool is_EM = proc_info.IsEM();
121 
122  if(is_CC && !is_delta) {
123  if((is_nu && is_p) || (is_nubar && is_n)) return 0;
124  }
125 
126  // Get baryon resonance parameters
127  int IR = utils::res::ResonanceIndex (resonance);
128  int LR = utils::res::OrbitalAngularMom (resonance);
129  double MR = utils::res::Mass (resonance);
130  double WR = utils::res::Width (resonance);
131  double NR = fNormBW?utils::res::BWNorm (resonance,fN0ResMaxNWidths,fN2ResMaxNWidths,fGnResMaxNWidths):1;
132 
133  // Following NeuGEN, avoid problems with underlying unphysical
134  // model assumptions by restricting the allowed W phase space
135  // around the resonance peak
136  if (fNormBW) {
137  if (W > MR + fN0ResMaxNWidths * WR && IR==0) return 0.;
138  else if (W > MR + fN2ResMaxNWidths * WR && IR==2) return 0.;
139  else if (W > MR + fGnResMaxNWidths * WR) return 0.;
140  }
141 
142  // Compute auxiliary & kinematical factors
143  double E = init_state.ProbeE(kRfHitNucRest);
144  double Mnuc = target.HitNucMass();
145  double W2 = TMath::Power(W, 2);
146  double Mnuc2 = TMath::Power(Mnuc, 2);
147  double k = 0.5 * (W2 - Mnuc2)/Mnuc;
148  double v = k - 0.5 * q2/Mnuc;
149  double v2 = TMath::Power(v, 2);
150  double Q2 = v2 - q2;
151  double Q = TMath::Sqrt(Q2);
152  double Eprime = E - v;
153  double U = 0.5 * (E + Eprime + Q) / E;
154  double V = 0.5 * (E + Eprime - Q) / E;
155  double U2 = TMath::Power(U, 2);
156  double V2 = TMath::Power(V, 2);
157  double UV = U*V;
158 
159 
160  //JN parameter from the KUZMIN et al.
161 
162  // bool is_RS = true;
163  bool is_KLN = false;
164  if(fKLN && is_CC) is_KLN=true;
165 
166  bool is_BRS = false;
167  if(fBRS && is_CC) is_BRS=true;
168 
169  double ml = interaction->FSPrimLepton()->Mass();
170  double Pl = TMath::Sqrt(Eprime*Eprime - ml*ml);
171 
172  double vstar = (Mnuc*v + q2)/W; //missing W
173  double Qstar = TMath::Sqrt(-q2 + vstar*vstar);
174  double sqrtq2 = TMath::Sqrt(-q2);
175  double a = 1. + 0.5*(W2-q2+Mnuc2)/Mnuc/W;
176 
177  double KNL_Alambda_plus = 0;
178  double KNL_Alambda_minus = 0;
179  double KNL_j0_plus = 0;
180  double KNL_j0_minus = 0;
181  double KNL_jx_plus = 0;
182  double KNL_jx_minus = 0;
183  double KNL_jy_plus = 0;
184  double KNL_jy_minus = 0;
185  double KNL_jz_plus = 0;
186  double KNL_jz_minus = 0;
187  double KNL_Qstar_plus =0;
188  double KNL_Qstar_minus =0;
189 
190  double KNL_K = Q/E/TMath::Sqrt(2*(-q2));
191 
192  double KNL_cL_plus = 0;
193  double KNL_cL_minus = 0;
194 
195  double KNL_cR_plus = 0;
196  double KNL_cR_minus = 0;
197 
198  double KNL_cS_plus = 0;
199  double KNL_cS_minus = 0;
200 
201  double KNL_vstar_plus = 0;
202  double KNL_vstar_minus = 0;
203 
204  if(is_CC && (is_KLN || is_BRS)){
205 
206  LOG("BSKLNBaseRESPXSec2014",pINFO) "costh1="<<costh;
207  costh = (q2 - ml*ml + 2.*E*Eprime)/2./E/Pl;
208  //ml=0;
209  LOG("BSKLNBaseRESPXSec2014",pINFO) "q2="<<q2<< "m2="<<ml*ml<<" 2.*E*Eprime="<<2.*E*Eprime<<" nom="<< (q2 - ml*ml + 2.*E*Eprime)<<" den="<<2.*E*Pl;
210  LOG("BSKLNBaseRESPXSec2014",pINFO) "costh2="<<costh;
211 
212  KNL_Alambda_plus = TMath::Sqrt(E*(Eprime - Pl));
213  KNL_Alambda_minus = TMath::Sqrt(E*(Eprime + Pl));
214  LOG("BSKLNBaseRESPXSec2014",pINFO)
215  << "\n+++++++++++++++++++++++ \n"
216  << "E="<<E << " K= "<<KNL_K << "\n"
217  << "El="<<Eprime<<" Pl="<<Pl<<" ml="<<ml << "\n"
218  << "W="<<W<<" Q="<<Q<<" q2="<<q2 << "\n"
219  << "A-="<<KNL_Alambda_minus<<" A+="<<KNL_Alambda_plus << "\n"
220  << "xxxxxxxxxxxxxxxxxxxxxxx";
221 
222  KNL_j0_plus = KNL_Alambda_plus /W * TMath::Sqrt(1 - costh) * (Mnuc - Eprime - Pl);
223  KNL_j0_minus = KNL_Alambda_minus/W * TMath::Sqrt(1 + costh) * (Mnuc - Eprime + Pl);
224 
225  KNL_jx_plus = KNL_Alambda_plus/ Q * TMath::Sqrt(1 + costh) * (Pl - E);
226  KNL_jx_minus = KNL_Alambda_minus/Q * TMath::Sqrt(1 - costh) * (Pl + E);
227 
228  KNL_jy_plus = KNL_Alambda_plus * TMath::Sqrt(1 + costh);
229  KNL_jy_minus = -KNL_Alambda_minus * TMath::Sqrt(1 - costh);
230 
231  KNL_jz_plus = KNL_Alambda_plus /W/Q * TMath::Sqrt(1 - costh) * ( (E + Pl)*(Mnuc -Eprime) + Pl*( E + 2*E*costh -Pl) );
232  KNL_jz_minus = KNL_Alambda_minus/W/Q * TMath::Sqrt(1 + costh) * ( (E - Pl)*(Mnuc -Eprime) + Pl*( -E + 2*E*costh -Pl) );
233 
234  if (is_nu || is_lminus) {
235  KNL_Qstar_plus = sqrtq2 * KNL_j0_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
236  KNL_Qstar_minus = sqrtq2 * KNL_j0_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
237  }
238 
239  else if (is_nubar || is_lplus){
240  KNL_Qstar_plus = sqrtq2 * KNL_j0_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
241  KNL_Qstar_minus = sqrtq2 * KNL_j0_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
242  }
243 
244  if (is_nu || is_lminus) {
245  KNL_vstar_plus = sqrtq2 * KNL_jz_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
246  KNL_vstar_minus = sqrtq2 * KNL_jz_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
247  }
248  else if (is_nubar || is_lplus) {
249  KNL_vstar_minus = sqrtq2 * KNL_jz_plus / TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
250  KNL_vstar_plus = sqrtq2 * KNL_jz_minus / TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
251  }
252 
253  if(is_nu || is_lminus){
254  KNL_cL_plus = TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus - KNL_jy_plus);
255  KNL_cL_minus = TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus - KNL_jy_minus);
256 
257  KNL_cR_plus = TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus + KNL_jy_plus);
258  KNL_cR_minus = TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus + KNL_jy_minus);
259 
260  KNL_cS_plus = KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_plus *KNL_j0_plus - KNL_jz_plus *KNL_jz_plus ) );
261  KNL_cS_minus = KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
262  }
263 
264  if (is_nubar || is_lplus) {
265  KNL_cL_plus = 1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus + KNL_jy_minus);
266  KNL_cL_minus = -1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus + KNL_jy_plus);
267 
268  KNL_cR_plus = 1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_minus - KNL_jy_minus);
269  KNL_cR_minus = -1 * TMath::Sqrt(0.5)* KNL_K * (KNL_jx_plus - KNL_jy_plus);
270 
271  KNL_cS_plus = -1 * KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_minus*KNL_j0_minus - KNL_jz_minus*KNL_jz_minus) );
272  KNL_cS_minus = 1 * KNL_K * TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
273  }
274  }
275 
276  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"j0-="<<KNL_j0_minus<<" j0+="<<KNL_j0_plus;
277  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"jx-="<<KNL_jx_minus<<" jx+="<<KNL_jx_plus;
278  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"jy-="<<KNL_jy_minus<<" jy+="<<KNL_jy_plus;
279  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"jz-="<<KNL_jz_minus<<" jz+="<<KNL_jz_plus;
280 
281  LOG("BSKLNBaseRESPXSec2014",pINFO) "sqrt2="<<sqrtq2<<" jz+=:"<<KNL_jz_plus<<" j0+="<<KNL_j0_plus<<" denom="<<TMath::Sqrt(TMath::Abs(KNL_j0_plus*KNL_j0_plus - KNL_jz_plus*KNL_jz_plus) );
282 
283  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"vstar-="<<KNL_vstar_minus<<" vstar+="<<KNL_vstar_plus;
284  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"Qstar-="<<KNL_Qstar_minus<<" Qstar+="<<KNL_Qstar_plus;
285 
286 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
287  LOG("BSKLNBaseRESPXSec2014", pDEBUG)
288  << "Kinematical params V = " << V << ", U = " << U;
289 #endif
290 
291  // Calculate RES Correction factor for neutrinos [ F.Ravndal, Nuovo Cim. A 18, 385 (1973) ]
292  double Go = TMath::Power( 1 - 0.25 * q2 / Mnuc2, 0.5 - IR ) ;
293 
294  // For EM, the correction factor is different [ F. Ravndal, Phys. Rev. D 4, 1466 (1971) ]
295  if( is_EM ) Go = TMath::Power( 1 - 0.25 * q2 / W2, 0.5*(1-IR) ) ;
296 
297  double GV = Go * TMath::Power( 1./(1-q2/fMv2), 2);
298  double GA = Go * TMath::Power( 1./(1-q2/fMa2), 2);
299 
300  if(fGVMiniBooNE){
301 
302  LOG("BSKLNBaseRESPXSec2014",pDEBUG) <<"Using new GV tuned to ANL and BNL data";
303  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"fCv3= " << fCv3 << ", fCv4= " << fCv4 << ", fCv51= " <<fCv51 << ", fCv52= " << fCv52;
304  double CV0 = 1./(1-q2/fMv2/4.);
305  double CV3 = fCv3 * CV0 * TMath::Power( 1-q2/fMv2,-2);
306  double CV4 = -1. * fCv4 * CV0 * TMath::Power( 1-q2/fMv2,-2);
307  double CV5 = fCv51* CV0 * TMath::Power( 1-q2/fMv2/fCv52, -2);
308 
309  double GV3 = 0.5 / TMath::Sqrt(3) * ( CV3 * (W + Mnuc)/Mnuc
310  + CV4 * (W2 + q2 -Mnuc2)/2./Mnuc2
311  + CV5 * (W2 - q2 -Mnuc2)/2./Mnuc2 );
312 
313  double GV1 = - 0.5 / TMath::Sqrt(3) * ( CV3 * (Mnuc2 -q2 +Mnuc*W)/W/Mnuc
314  + CV4 * (W2 +q2 - Mnuc2)/2./Mnuc2
315  + CV5 * (W2 -q2 - Mnuc2)/2./Mnuc2 );
316 
317  GV = 0.5 * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), 0.5-IR)
318  * TMath::Sqrt( 3 * GV3*GV3 + GV1*GV1);
319 
320  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"GV= " <<GV << " CV3= " <<CV3 << " CV4= " << CV4 << " CV5= " << CV5 << " GV3= " << GV3 << " GV1= " <<GV1;
321  } else if(fGVSaritaSchwinger){
322  // PhysRevD.77.053001
323  LOG("BSKLNBaseRESPXSec2014",pDEBUG) <<"Using GV Sarita-Schwinger model";
324  double CV0 = 1./(1-q2/fMv2/4.);
325  double CV3 = fCv3 * CV0 * TMath::Power( 1-q2/fMv2,-2);
326  double CV4 = -1. * fCv4 * CV0 * TMath::Power( 1-q2/fMv2,-2);
327  double CV5 = fCv51 * CV0 * TMath::Power( 1-q2/fMv2/fCv52, -2);
328 
329  double GV3 = 0.5 / TMath::Sqrt(3) * ( CV3 * (W + Mnuc)/Mnuc
330  + CV4 * (W2 + q2 -Mnuc2)/2./Mnuc2
331  + CV5 * (W2 - q2 -Mnuc2)/2./Mnuc2 );
332 
333  double GV1 = - 0.5 / TMath::Sqrt(3) * ( CV3 * (Mnuc2 -q2 +Mnuc*W)/W/Mnuc
334  + CV4 * (W2 +q2 - Mnuc2)/2./Mnuc2
335  + CV5 * (W2 -q2 - Mnuc2)/2./Mnuc2 );
336 
337  GV = 0.5 * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + Mnuc), -IR);
338 
339  if( is_EM ) GV = 0.5 * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), -0.5*IR);
340  GV *= TMath::Sqrt( 3 * GV3*GV3 + GV1*GV1);
341 
342  } else {
343  LOG("BSKLNBaseRESPXSec2014",pDEBUG << "Using dipole parametrization for GV") ;
344  }
345 
346  if(fGAMiniBooNE){
347  LOG("BSKLNBaseRESPXSec2014",pDEBUG) << "Using new GA tuned to ANL and BNL data";
348 
349  double CA5 = fCa50 * TMath::Power( 1./(1-q2/fMa2), 2);
350  // GA = 0.5 * TMath::Sqrt(3.) * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), 0.5-IR) * (1- (W2 +q2 -Mnuc2)/8./Mnuc2) * CA5/fZeta;
351  GA = 0.5 * TMath::Sqrt(3.) * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), 0.5-IR) * (1- (W2 +q2 -Mnuc2)/8./Mnuc2) * CA5;
352  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"GA= " <<GA << " CA50= " <<fCa50 << " C5A= " <<CA5;
353 
354  } else if(fGASaritaSchwinger){
355  LOG("BSKLNBaseRESPXSec2014",pDEBUG) << "Using GA Rarita-Schwinger model";
356 
357  double CA5 = fCa50 * TMath::Power( 1./(1-q2/fMa2), 2) * ( 1./(1 - fcII * q2/fMb2) );
358  GA = 0.5 * TMath::Sqrt(3.) * TMath::Power( 1 - q2/(Mnuc + W)/(Mnuc + W), 0.5) * TMath::Power( 1 - q2/(4*Mnuc2), -IR) * (1- (W2 +q2 -Mnuc2)/8./Mnuc2) * CA5;
359 
360  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"GA= " <<GA << " C5A= " <<CA5;
361 
362  } else {
363  LOG("BSKLNBaseRESPXSec2014",pDEBUG << "Using dipole parametrization for GA") ;
364  }
365 
366  if(is_EM) {
367  GA = 0.; // zero the axial term for EM scattering
368  }
369 
370  double d = TMath::Power(W+Mnuc,2.) - q2;
371  double sq2omg = TMath::Sqrt(2./fOmega);
372  double nomg = IR * fOmega;
373  double mq_w = Mnuc*Q/W;
374 
375  fFKR.Lamda = sq2omg * mq_w;
376  fFKR.Tv = GV / (3.*W*sq2omg);
377  fFKR.Rv = kSqrt2 * mq_w*(W+Mnuc)*GV / d;
378  fFKR.S = (-q2/Q2) * (3*W*Mnuc + q2 - Mnuc2) * GV / (6*Mnuc2);
379  fFKR.Ta = (2./3.) * (fZeta/sq2omg) * mq_w * GA / d;
380  fFKR.Ra = (kSqrt2/6.) * fZeta * (GA/W) * (W+Mnuc + 2*nomg*W/d );
381  fFKR.B = fZeta/(3.*W*sq2omg) * (1 + (W2-Mnuc2+q2)/ d) * GA;
382  fFKR.C = fZeta/(6.*Q) * (W2 - Mnuc2 + nomg*(W2-Mnuc2+q2)/d) * (GA/Mnuc);
383  fFKR.R = fFKR.Rv;
384  fFKR.Rplus = - (fFKR.Rv + fFKR.Ra);
385  fFKR.Rminus = - (fFKR.Rv - fFKR.Ra);
386  fFKR.T = fFKR.Tv;
387  fFKR.Tplus = - (fFKR.Tv + fFKR.Ta);
388  fFKR.Tminus = - (fFKR.Tv - fFKR.Ta);
389 
390  //JN KNL
391  double KNL_S_plus = 0;
392  double KNL_S_minus = 0;
393  double KNL_B_plus = 0;
394  double KNL_B_minus = 0;
395  double KNL_C_plus = 0;
396  double KNL_C_minus = 0;
397 
398  if(is_CC && is_KLN){
399  KNL_S_plus = (KNL_vstar_plus*vstar - KNL_Qstar_plus *Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2; //possibly missing minus sign ()
400  KNL_S_minus = (KNL_vstar_minus*vstar - KNL_Qstar_minus*Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2;
401 
402  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"KNL S= " <<KNL_S_plus<<"\t"<<KNL_S_minus<<"\t"<<fFKR.S;
403 
404  KNL_B_plus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_plus + KNL_vstar_plus *Qstar/a/Mnuc ) * GA;
405  KNL_B_minus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_minus + KNL_vstar_minus*Qstar/a/Mnuc ) * GA;
406  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"KNL B= " <<KNL_B_plus<<"\t"<<KNL_B_minus<<"\t"<<fFKR.B;
407 
408  KNL_C_plus = ( (KNL_Qstar_plus*Qstar - KNL_vstar_plus*vstar ) * ( 1./3. + vstar/a/Mnuc)
409  + KNL_vstar_plus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
410 
411  KNL_C_minus = ( (KNL_Qstar_minus*Qstar - KNL_vstar_minus*vstar ) * ( 1./3. + vstar/a/Mnuc)
412  + KNL_vstar_minus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
413 
414  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"KNL C= "<<KNL_C_plus<<"\t"<<KNL_C_minus<<"\t"<<fFKR.C;
415  }
416  double BRS_S_plus = 0;
417  double BRS_S_minus = 0;
418  double BRS_B_plus = 0;
419  double BRS_B_minus = 0;
420  double BRS_C_plus = 0;
421  double BRS_C_minus = 0;
422 
423 
424  if(is_CC && is_BRS){
425 
426  KNL_S_plus = (KNL_vstar_plus*vstar - KNL_Qstar_plus *Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2;
427  KNL_S_minus = (KNL_vstar_minus*vstar - KNL_Qstar_minus*Qstar )* (Mnuc2 -q2 - 3*W*Mnuc ) * GV / (6*Mnuc2)/Q2;
428 
429 
430  KNL_B_plus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_plus + KNL_vstar_plus *Qstar/a/Mnuc ) * GA;
431  KNL_B_minus = fZeta/(3.*W*sq2omg)/Qstar * (KNL_Qstar_minus + KNL_vstar_minus*Qstar/a/Mnuc ) * GA;
432 
433 
434  KNL_C_plus = ( (KNL_Qstar_plus*Qstar - KNL_vstar_plus*vstar ) * ( 1./3. + vstar/a/Mnuc)
435  + KNL_vstar_plus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
436 
437  KNL_C_minus = ( (KNL_Qstar_minus*Qstar - KNL_vstar_minus*vstar ) * ( 1./3. + vstar/a/Mnuc)
438  + KNL_vstar_minus*(2./3.*W +q2/a/Mnuc + nomg/3./a/Mnuc) )* fZeta * (GA/2./W/Qstar);
439 
440  BRS_S_plus = KNL_S_plus;
441  BRS_S_minus = KNL_S_minus;
442  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"BRS S= " <<KNL_S_plus<<"\t"<<KNL_S_minus<<"\t"<<fFKR.S;
443 
444  BRS_B_plus = KNL_B_plus + fZeta*GA/2./W/Qstar*( KNL_Qstar_plus*vstar - KNL_vstar_plus*Qstar)
445  *( 2./3 /sq2omg *(vstar + Qstar*Qstar/Mnuc/a))/(kPionMass2 -q2);
446 
447  BRS_B_minus = KNL_B_minus + fZeta*GA/2./W/Qstar*( KNL_Qstar_minus*vstar - KNL_vstar_minus*Qstar)
448  *( 2./3 /sq2omg *(vstar + Qstar*Qstar/Mnuc/a))/(kPionMass2 -q2);
449  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"BRS B= " <<KNL_B_plus<<"\t"<<KNL_B_minus<<"\t"<<fFKR.B;
450 
451  BRS_C_plus = KNL_C_plus + fZeta*GA/2./W/Qstar*( KNL_Qstar_plus*vstar - KNL_vstar_plus*Qstar)
452  * Qstar*(2./3.*W +q2/Mnuc/a +nomg/3./a/Mnuc)/(kPionMass2 -q2);
453 
454  BRS_C_minus = KNL_C_minus + fZeta*GA/2./W/Qstar*( KNL_Qstar_minus*vstar - KNL_vstar_minus*Qstar)
455  * Qstar*(2./3.*W +q2/Mnuc/a +nomg/3./a/Mnuc)/(kPionMass2 -q2);
456  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"BRS C= " <<KNL_C_plus<<"\t"<<KNL_C_minus<<"\t"<<fFKR.C;
457  }
458 
459 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
460  LOG("FKR", pDEBUG)
461  << "FKR params for RES = " << resname << " : " << fFKR;
462 #endif
463 
464  // Calculate the Rein-Sehgal Helicity Amplitudes
465  double sigL_minus = 0;
466  double sigR_minus = 0;
467  double sigS_minus = 0;
468 
469  double sigL_plus = 0;
470  double sigR_plus = 0;
471  double sigS_plus = 0;
472 
473  const RSHelicityAmplModelI * hamplmod = 0;
474  const RSHelicityAmplModelI * hamplmod_KNL_minus = 0;
475  const RSHelicityAmplModelI * hamplmod_KNL_plus = 0;
476  const RSHelicityAmplModelI * hamplmod_BRS_minus = 0;
477  const RSHelicityAmplModelI * hamplmod_BRS_plus = 0;
478 
479 
480  double g2 = kGF2; // NC
481  if(is_CC) g2 *= fVud2;
482 
483  if(is_EM)
484  {
485  double q4 = q2*q2;
486  g2 = 8*kAem2*kPi2/q4;
487  }
488 
489  double sig0 = 0.125*(g2/kPi)*(-q2/Q2)*(W/Mnuc);
490  double scLR = W/Mnuc;
491  double scS = (Mnuc/W)*(-Q2/q2);
492 
493  double sigL =0;
494  double sigR =0;
495  double sigS =0;
496 
497  double sigRSL =0;
498  double sigRSR =0;
499  double sigRSS =0;
500 
501  if(is_CC && !(is_KLN || is_BRS) ) {
502 
503  hamplmod = fHAmplModelCC;
504  }
505  else
506  if(is_NC) {
507  if (is_p) { hamplmod = fHAmplModelNCp;}
508  else { hamplmod = fHAmplModelNCn;}
509  }
510  else
511  if(is_EM) {
512  if (is_p) { hamplmod = fHAmplModelEMp;}
513  else { hamplmod = fHAmplModelEMn;}
514  }
515  else
516  if(is_CC && is_KLN ){
517  fFKR.S = KNL_S_minus; //2 times fFKR.S?
518  fFKR.B = KNL_B_minus;
519  fFKR.C = KNL_C_minus;
520 
521  hamplmod_KNL_minus = fHAmplModelCC;
522 
523  assert(hamplmod_KNL_minus);
524 
525  const RSHelicityAmpl & hampl_KNL_minus = hamplmod_KNL_minus->Compute(resonance, fFKR);
526 
527  sigL_minus = (hampl_KNL_minus.Amp2Plus3 () + hampl_KNL_minus.Amp2Plus1 ());
528  sigR_minus = (hampl_KNL_minus.Amp2Minus3() + hampl_KNL_minus.Amp2Minus1());
529  sigS_minus = (hampl_KNL_minus.Amp20Plus () + hampl_KNL_minus.Amp20Minus());
530 
531 
532  fFKR.S = KNL_S_plus;
533  fFKR.B = KNL_B_plus;
534  fFKR.C = KNL_C_plus;
535  hamplmod_KNL_plus = fHAmplModelCC;
536  assert(hamplmod_KNL_plus);
537 
538  const RSHelicityAmpl & hampl_KNL_plus = hamplmod_KNL_plus->Compute(resonance, fFKR);
539 
540  sigL_plus = (hampl_KNL_plus.Amp2Plus3 () + hampl_KNL_plus.Amp2Plus1 ());
541  sigR_plus = (hampl_KNL_plus.Amp2Minus3() + hampl_KNL_plus.Amp2Minus1());
542  sigS_plus = (hampl_KNL_plus.Amp20Plus () + hampl_KNL_plus.Amp20Minus());
543 
544  }
545  else
546  if(is_CC && is_BRS ){
547  fFKR.S = BRS_S_minus;
548  fFKR.B = BRS_B_minus;
549  fFKR.C = BRS_C_minus;
550 
551  hamplmod_BRS_minus = fHAmplModelCC;
552  assert(hamplmod_BRS_minus);
553 
554  const RSHelicityAmpl & hampl_BRS_minus = hamplmod_BRS_minus->Compute(resonance, fFKR);
555 
556  sigL_minus = (hampl_BRS_minus.Amp2Plus3 () + hampl_BRS_minus.Amp2Plus1 ());
557  sigR_minus = (hampl_BRS_minus.Amp2Minus3() + hampl_BRS_minus.Amp2Minus1());
558  sigS_minus = (hampl_BRS_minus.Amp20Plus () + hampl_BRS_minus.Amp20Minus());
559 
560  fFKR.S = BRS_S_plus;
561  fFKR.B = BRS_B_plus;
562  fFKR.C = BRS_C_plus;
563  hamplmod_BRS_plus = fHAmplModelCC;
564  assert(hamplmod_BRS_plus);
565 
566  const RSHelicityAmpl & hampl_BRS_plus = hamplmod_BRS_plus->Compute(resonance, fFKR);
567 
568  sigL_plus = (hampl_BRS_plus.Amp2Plus3 () + hampl_BRS_plus.Amp2Plus1 ());
569  sigR_plus = (hampl_BRS_plus.Amp2Minus3() + hampl_BRS_plus.Amp2Minus1());
570  sigS_plus = (hampl_BRS_plus.Amp20Plus () + hampl_BRS_plus.Amp20Minus());
571  }
572 
573  // Compute the cross section
574  if(is_KLN || is_BRS) {
575 
576  sigL_minus *= scLR;
577  sigR_minus *= scLR;
578  sigS_minus *= scS;
579  sigL_plus *= scLR;
580  sigR_plus *= scLR;
581  sigS_plus *= scS;
582 
583  LOG("BSKLNBaseRESPXSec2014", pINFO)
584  << "sL,R,S minus = " << sigL_minus << "," << sigR_minus << "," << sigS_minus;
585  LOG("BSKLNBaseRESPXSec2014", pINFO)
586  << "sL,R,S plus = " << sigL_plus << "," << sigR_plus << "," << sigS_plus;
587  }
588  else {
589  assert(hamplmod);
590 
591  const RSHelicityAmpl & hampl = hamplmod->Compute(resonance, fFKR);
592 
593  sigL = scLR* (hampl.Amp2Plus3 () + hampl.Amp2Plus1 ());
594  sigR = scLR* (hampl.Amp2Minus3() + hampl.Amp2Minus1());
595  sigS = scS * (hampl.Amp20Plus () + hampl.Amp20Minus());
596  }
597 
598 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
599  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{0} = " << sig0;
600  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{L} = " << sigL;
601  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{R} = " << sigR;
602  LOG("BSKLNBaseRESPXSec2014", pDEBUG) << "sig_{S} = " << sigS;
603 #endif
604 
605  double xsec = 0.0;
606 
607  if(is_KLN || is_BRS) {
608  xsec = TMath::Power(KNL_cL_minus,2)*sigL_minus + TMath::Power(KNL_cL_plus,2)*sigL_plus
609  + TMath::Power(KNL_cR_minus,2)*sigR_minus + TMath::Power(KNL_cR_plus,2)*sigR_plus
610  + TMath::Power(KNL_cS_minus,2)*sigS_minus + TMath::Power(KNL_cS_plus,2)*sigS_plus;
611  xsec *=sig0;
612 
613  LOG("BSKLNBaseRESPXSec2014",pINFO) << "A-="<<KNL_Alambda_minus<<" A+="<<KNL_Alambda_plus;
614  // protect against sigRSR=sigRSL=sigRSS=0
615  LOG("BSKLNBaseRESPXSec2014",pINFO) <<q2<<"\t"<<xsec<<"\t"<<sig0*(V2*sigR + U2*sigL + 2*UV*sigS)<<"\t"<<xsec/TMath::Max(sig0*(V2*sigRSR + U2*sigRSL + 2*UV*sigRSS),1.0e-100);
616  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"fFKR.B="<<fFKR.B<<" fFKR.C="<<fFKR.C<<" fFKR.S="<<fFKR.S;
617  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"CL-="<<TMath::Power(KNL_cL_minus,2)<<" CL+="<<TMath::Power(KNL_cL_plus,2)<<" U2="<<U2;
618  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"SL-="<<sigL_minus<<" SL+="<<sigL_plus<<" SL="<<sigRSL;
619 
620  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"CR-="<<TMath::Power(KNL_cR_minus,2)<<" CR+="<<TMath::Power(KNL_cR_plus,2)<<" V2="<<V2;
621  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"SR-="<<sigR_minus<<" SR+="<<sigR_plus<<" sR="<<sigRSR;
622 
623  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"CS-="<<TMath::Power(KNL_cS_minus,2)<<" CS+="<<TMath::Power(KNL_cS_plus,2)<<" UV="<<UV;
624  LOG("BSKLNBaseRESPXSec2014",pINFO) <<"SS-="<<sigL_minus<<" SS+="<<sigS_plus<<" sS="<<sigRSS;
625  }
626  else {
627  if (is_nu || is_lminus) {
628  xsec = sig0*(V2*sigR + U2*sigL + 2*UV*sigS);
629  }
630  else
631  if (is_nubar || is_lplus) {
632  xsec = sig0*(U2*sigR + V2*sigL + 2*UV*sigS);
633  }
634  xsec = TMath::Max(0.,xsec);
635  }
636  double mult = 1.0;
637  if ( is_CC && is_delta ) {
638  if ( (is_nu && is_p) || (is_nubar && is_n) ) mult=3.0;
639  }
640  xsec *= mult;
641 
642  // Check whether the cross section is to be weighted with a Breit-Wigner distribution
643  // (default: true)
644  double bw = 1.0;
645  if ( fWghtBW ) {
646  bw = utils::bwfunc::BreitWignerL(W,LR,MR,WR,NR);
647  }
648 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
649  LOG("BSKLNBaseRESPXSec2014", pDEBUG)
650  << "BreitWigner(RES=" << resname << ", W=" << W << ") = " << bw;
651 #endif
652  xsec *= bw;
653 
654 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
655  LOG("BSKLNBaseRESPXSec2014", pINFO)
656  << "\n d2xsec/dQ2dW" << "[" << interaction->AsString()
657  << "](W=" << W << ", q2=" << q2 << ", E=" << E << ") = " << xsec;
658 #endif
659 
660  // The algorithm computes d^2xsec/dWdQ2
661  // Check whether variable tranformation is needed
662  if ( kps != kPSWQ2fE ) {
663  double J = utils::kinematics::Jacobian(interaction,kPSWQ2fE,kps);
664  xsec *= J;
665  }
666 
667  // Apply given scaling factor
668  if (is_CC) { xsec *= fXSecScaleCC; }
669  else if (is_NC) { xsec *= fXSecScaleNC; }
670  else if (is_EM) { xsec *= fXSecScaleEM; }
671 
672  // If requested return the free nucleon xsec even for input nuclear tgt
673  if ( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
674 
675  int Z = target.Z();
676  int A = target.A();
677  int N = A-Z;
678 
679  // Take into account the number of scattering centers in the target
680  int NNucl = (is_p) ? Z : N;
681  xsec*=NNucl; // nuclear xsec (no nuclear suppression factor)
682 
683  if ( fUsePauliBlocking && A!=1 )
684  {
685  // Calculation of Pauli blocking according references:
686  //
687  // [1] S.L. Adler, S. Nussinov, and E.A. Paschos, "Nuclear
688  // charge exchange corrections to leptonic pion production
689  // in the (3,3) resonance region," Phys. Rev. D 9 (1974)
690  // 2125-2143 [Erratum Phys. Rev. D 10 (1974) 1669].
691  // [2] J.Y. Yu, "Neutrino interactions and nuclear effects in
692  // oscillation experiments and the nonperturbative disper-
693  // sive sector in strong (quasi-)abelian fields," Ph. D.
694  // Thesis, Dortmund U., Dortmund, 2002 (unpublished).
695  // [3] E.A. Paschos, J.Y. Yu, and M. Sakuda, "Neutrino pro-
696  // duction of resonances," Phys. Rev. D 69 (2004) 014013
697  // [arXiv: hep-ph/0308130].
698 
699  double P_Fermi = 0.0;
700 
701  // Maximum value of Fermi momentum of target nucleon (GeV)
702  if ( A<6 || ! fUseRFGParametrization )
703  {
704  // look up the Fermi momentum for this target
705  FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance();
706  const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
707  P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucpdgc);
708  }
709  else {
710  // define the Fermi momentum for this target
711  P_Fermi = utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(target);
712  // correct the Fermi momentum for the struck nucleon
713  if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); }
714  else { P_Fermi *= TMath::Power( 2.*N/A, 1./3); }
715  }
716 
717  double FactorPauli_RES = 1.0;
718 
719  double k0 = 0., q = 0., q0 = 0.;
720 
721  if (P_Fermi > 0.)
722  {
723  k0 = (W2-Mnuc2-Q2)/(2*W);
724  k = TMath::Sqrt(k0*k0+Q2); // previous value of k is overridden
725  q0 = (W2-Mnuc2+kPionMass2)/(2*W);
726  q = TMath::Sqrt(q0*q0-kPionMass2);
727 
728  if ( 2*P_Fermi < k-q )
729  FactorPauli_RES = 1.0;
730  if ( 2*P_Fermi >= k+q )
731  FactorPauli_RES = ((3*k*k+q*q)/(2*P_Fermi)-(5*TMath::Power(k,4)+TMath::Power(q,4)+10*k*k*q*q)/(40*TMath::Power(P_Fermi,3)))/(2*k);
732  if ( 2*P_Fermi >= k-q && 2*P_Fermi <= k+q )
733  FactorPauli_RES = ((q+k)*(q+k)-4*P_Fermi*P_Fermi/5-TMath::Power(k-q, 3)/(2*P_Fermi)+TMath::Power(k-q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*q*k);
734  }
735 
736  xsec *= FactorPauli_RES;
737  }
738  return xsec;
739 }
740 //____________________________________________________________________________
741 double BSKLNBaseRESPXSec2014::Integral(const Interaction * interaction) const
742 {
743  double xsec = fXSecIntegrator->Integrate(this,interaction);
744  return xsec;
745 }
746 //____________________________________________________________________________
747 bool BSKLNBaseRESPXSec2014::ValidProcess(const Interaction * interaction) const
748 {
749  if(interaction->TestBit(kISkipProcessChk)) return true;
750 
751  const InitialState & init_state = interaction->InitState();
752  const ProcessInfo & proc_info = interaction->ProcInfo();
753  const XclsTag & xcls = interaction->ExclTag();
754 
755  if(!proc_info.IsResonant()) return false;
756  if(!xcls.KnownResonance()) return false;
757 
758  int hitnuc = init_state.Tgt().HitNucPdg();
759  bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc));
760 
761  if (!is_pn) return false;
762 
763  int probe = init_state.ProbePdg();
764  bool is_weak = proc_info.IsWeak();
765  bool is_em = proc_info.IsEM();
766  bool nu_weak = (pdg::IsNeutralLepton(probe) && is_weak);
767  bool l_em = (pdg::IsChargedLepton(probe) && is_em );
768 
769  if (!nu_weak && !l_em) return false;
770 
771  return true;
772 }
773 //____________________________________________________________________________
774 void BSKLNBaseRESPXSec2014::Configure(const Registry & config)
775 {
776  Algorithm::Configure(config);
777  this->LoadConfig();
778 }
779 //____________________________________________________________________________
780 void BSKLNBaseRESPXSec2014::Configure(string config)
781 {
782  Algorithm::Configure(config);
783  this->LoadConfig();
784 }
785 //____________________________________________________________________________
786 void BSKLNBaseRESPXSec2014::LoadConfig(void)
787 {
788  // Cross section scaling factors
789  this->GetParam( "RES-CC-XSecScale", fXSecScaleCC ) ;
790  this->GetParam( "RES-NC-XSecScale", fXSecScaleNC ) ;
791  this->GetParam( "RES-EM-XSecScale", fXSecScaleEM ) ;
792 
793  // Load all configuration data or set defaults
794 
795  this->GetParam( "RES-Zeta" , fZeta ) ;
796  this->GetParam( "RES-Omega" , fOmega ) ;
797  this->GetParam( "minibooneGA", fGAMiniBooNE ) ;
798  this->GetParam( "minibooneGV", fGVMiniBooNE ) ;
799  this->GetParam( "GASaritaSchwinger", fGASaritaSchwinger ) ;
800  this->GetParam( "GVSaritaSchwinger", fGVSaritaSchwinger ) ;
801 
802  double fermi_constant ;
803  this->GetParam( "FermiConstant", fermi_constant ) ;
804  fFermiConstant2 = fermi_constant * fermi_constant ;
805 
806  double alpha ;
807  this->GetParam( "FineStructureConstant", alpha ) ;
808  fFineStructure2 = alpha * alpha ;
809 
810  double ma, mv ;
811  this->GetParam( "RES-Ma", ma ) ;
812  this->GetParam( "RES-Mv", mv ) ;
813  fMa2 = TMath::Power(ma,2);
814  fMv2 = TMath::Power(mv,2);
815 
816  // Additional parameters used for the Sarita-Schwinger parameterization of GV and GA
817  // PhysRevD.77.053001
818  double mb;
819  this->GetParam( "GVCAL-Cv3" , fCv3) ;
820  this->GetParam( "GVCAL-Cv4" , fCv4) ;
821  this->GetParam( "GVCAL-Cv51" , fCv51) ;
822  this->GetParam( "GVCAL-Cv52" , fCv52) ;
823  this->GetParam( "RES-CA50", fCa50 ) ;
824  this->GetParamDef( "GAcII", fcII, 0. ) ;
825  this->GetParamDef( "RES-Mb", mb, 1. ) ;
826  fMb2 = TMath::Power(mb,2);
827 
828  this->GetParamDef( "BreitWignerWeight", fWghtBW, true ) ;
829  this->GetParamDef( "BreitWignerNorm", fNormBW, true);
830  double Vud;
831  this->GetParam("CKM-Vud", Vud );
832  fVud2 = TMath::Power( Vud, 2 );
833  this->GetParam("FermiMomentumTable", fKFTable);
834  this->GetParam("RFG-UseParametrization", fUseRFGParametrization);
835  this->GetParam("UsePauliBlockingForRES", fUsePauliBlocking);
836 
837  // Load all the sub-algorithms needed
838 
839  fHAmplModelCC = 0;
840  fHAmplModelNCp = 0;
841  fHAmplModelNCn = 0;
842  fHAmplModelEMp = 0;
843  fHAmplModelEMn = 0;
844 
845  AlgFactory * algf = AlgFactory::Instance();
846 
847  fHAmplModelCC = dynamic_cast<const RSHelicityAmplModelI *> (
848  algf->GetAlgorithm("genie::RSHelicityAmplModelCC","Default"));
849  fHAmplModelNCp = dynamic_cast<const RSHelicityAmplModelI *> (
850  algf->GetAlgorithm("genie::RSHelicityAmplModelNCp","Default"));
851  fHAmplModelNCn = dynamic_cast<const RSHelicityAmplModelI *> (
852  algf->GetAlgorithm("genie::RSHelicityAmplModelNCn","Default"));
853  fHAmplModelEMp = dynamic_cast<const RSHelicityAmplModelI *> (
854  algf->GetAlgorithm("genie::RSHelicityAmplModelEMp","Default"));
855  fHAmplModelEMn = dynamic_cast<const RSHelicityAmplModelI *> (
856  algf->GetAlgorithm("genie::RSHelicityAmplModelEMn","Default"));
857 
858  assert( fHAmplModelCC );
859  assert( fHAmplModelNCp );
860  assert( fHAmplModelNCn );
861  assert( fHAmplModelEMp );
862  assert( fHAmplModelEMn );
863 
864  // Use algorithm within a DIS/RES join scheme. If yes get Wcut
865  this->GetParam( "UseDRJoinScheme", fUsingDisResJoin ) ;
866  fWcut = 999999;
867  if(fUsingDisResJoin) {
868  this->GetParam( "Wcut", fWcut ) ;
869  }
870 
871  // NeuGEN limits in the allowed resonance phase space:
872  // W < min{ Wmin(physical), (res mass) + x * (res width) }
873  // It limits the integration area around the peak and avoids the
874  // problem with huge xsec increase at low Q2 and high W.
875  // In correspondence with Hugh, Rein said that the underlying problem
876  // are unphysical assumptions in the model.
877  this->GetParamDef( "MaxNWidthForN2Res", fN2ResMaxNWidths, 2.0 ) ;
878  this->GetParamDef( "MaxNWidthForN0Res", fN0ResMaxNWidths, 6.0 ) ;
879  this->GetParamDef( "MaxNWidthForGNRes", fGnResMaxNWidths, 4.0 ) ;
880 
881  // Load the differential cross section integrator
882  fXSecIntegrator =
883  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
884  assert(fXSecIntegrator);
885 }
886 //____________________________________________________________________________
genie::utils::res::IsDelta
bool IsDelta(Resonance_t res)
is it a Delta resonance?
Definition: BaryonResUtils.cxx:398
genie::ProcessInfo::IsResonant
bool IsResonant(void) const
Definition: ProcessInfo.cxx:99
genie::BSKLNBaseRESPXSec2014::fNormBW
bool fNormBW
normalize resonance breit-wigner to 1?
Definition: BSKLNBaseRESPXSec2014.h:91
genie::XSecAlgorithmI
Cross Section Calculation Interface.
Definition: XSecAlgorithmI.h:27
genie::BSKLNBaseRESPXSec2014::fKFTable
string fKFTable
table of Fermi momentum (kF) constants for various nuclei
Definition: BSKLNBaseRESPXSec2014.h:103
genie::BSKLNBaseRESPXSec2014::fOmega
double fOmega
FKR parameter Omega.
Definition: BSKLNBaseRESPXSec2014.h:93
genie::Kinematics::W
double W(bool selected=false) const
Definition: Kinematics.cxx:157
genie::BSKLNBaseRESPXSec2014::Configure
void Configure(const Registry &config)
Definition: BSKLNBaseRESPXSec2014.cxx:774
genie::ProcessInfo::IsWeak
bool IsWeak(void) const
Definition: ProcessInfo.cxx:198
genie::ProcessInfo::IsWeakCC
bool IsWeakCC(void) const
Definition: ProcessInfo.cxx:203
genie::pdg::IsNeutrino
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
genie::BSKLNBaseRESPXSec2014::fUsingDisResJoin
bool fUsingDisResJoin
use a DIS/RES joining scheme?
Definition: BSKLNBaseRESPXSec2014.h:98
genie::BSKLNBaseRESPXSec2014::fXSecScaleNC
double fXSecScaleNC
external NC xsec scaling factor
Definition: BSKLNBaseRESPXSec2014.h:108
genie::utils::mec::J
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
genie::FKR::Rminus
double Rminus
Definition: FKR.h:50
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::FKR::Ra
double Ra
Definition: FKR.h:42
genie::RSHelicityAmpl::Amp2Plus3
double Amp2Plus3(void) const
Definition: RSHelicityAmpl.h:66
genie::Target::A
int A(void) const
Definition: Target.h:70
genie::RSHelicityAmpl::Amp2Minus3
double Amp2Minus3(void) const
Definition: RSHelicityAmpl.h:65
genie::RSHelicityAmplModelI::Compute
virtual const RSHelicityAmpl & Compute(Resonance_t res, const FKR &fkr) const =0
genie::XclsTag::KnownResonance
bool KnownResonance(void) const
Definition: XclsTag.h:68
genie::Target::HitNucMass
double HitNucMass(void) const
Definition: Target.cxx:233
genie::BSKLNBaseRESPXSec2014::fN0ResMaxNWidths
double fN0ResMaxNWidths
limits allowed phase space for n=0 res
Definition: BSKLNBaseRESPXSec2014.h:101
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::FKR::Lamda
double Lamda
Definition: FKR.h:37
genie::pdg::IsChargedLepton
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
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::RSHelicityAmpl::Amp2Plus1
double Amp2Plus1(void) const
Definition: RSHelicityAmpl.h:64
genie::RSHelicityAmpl::Amp2Minus1
double Amp2Minus1(void) const
return |helicity amplitude|^2
Definition: RSHelicityAmpl.h:63
genie::utils::bwfunc::BreitWignerL
double BreitWignerL(double W, int L, double mass, double width0, double norm)
Definition: BWFunc.cxx:99
genie::KinePhaseSpace_t
enum genie::EKinePhaseSpace KinePhaseSpace_t
genie::utils::res::BWNorm
double BWNorm(Resonance_t res, double N0ResMaxNWidths=6, double N2ResMaxNWidths=2, double GnResMaxNWidths=4)
breit-wigner normalization factor
Definition: BaryonResUtils.cxx:490
genie::Resonance_t
enum genie::EResonance Resonance_t
genie::BSKLNBaseRESPXSec2014::fHAmplModelEMp
const RSHelicityAmplModelI * fHAmplModelEMp
Definition: BSKLNBaseRESPXSec2014.h:83
genie::Interaction::AsString
string AsString(void) const
Definition: Interaction.cxx:246
genie::XclsTag
Contains minimal information for tagging exclusive processes.
Definition: XclsTag.h:39
genie::BSKLNBaseRESPXSec2014::fHAmplModelCC
const RSHelicityAmplModelI * fHAmplModelCC
Definition: BSKLNBaseRESPXSec2014.h:80
genie::utils::kinematics::W
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
genie::BSKLNBaseRESPXSec2014::fVud2
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
Definition: BSKLNBaseRESPXSec2014.h:97
genie::pdg::IsNeutron
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
genie::pdg::IsPosChargedLepton
bool IsPosChargedLepton(int pdgc)
Definition: PDGUtils.cxx:148
genie::Interaction
Summary information for an interaction.
Definition: Interaction.h:56
genie::FKR::Tv
double Tv
Definition: FKR.h:38
genie::XSecAlgorithmI::ValidKinematics
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
Definition: XSecAlgorithmI.cxx:40
genie::Kinematics::q2
double q2(bool selected=false) const
Definition: Kinematics.cxx:141
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::Kinematics::FSLeptonP4
const TLorentzVector & FSLeptonP4(void) const
Definition: Kinematics.h:65
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::BSKLNBaseRESPXSec2014::fWghtBW
bool fWghtBW
weight with resonance breit-wigner?
Definition: BSKLNBaseRESPXSec2014.h:90
genie::units::A
static constexpr double A
Definition: Units.h:74
genie::FermiMomentumTablePool::GetTable
const FermiMomentumTable * GetTable(string name)
Definition: FermiMomentumTablePool.cxx:76
genie::ProcessInfo
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
a
const double a
Definition: gUpMuFluxGen.cxx:164
genie::units::mb
static constexpr double mb
Definition: Units.h:79
genie::FKR::T
double T
Definition: FKR.h:46
genie::FKR::Rv
double Rv
Definition: FKR.h:39
genie::pdg::IsAntiNeutrino
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
genie::BSKLNBaseRESPXSec2014::fXSecScaleCC
double fXSecScaleCC
external CC xsec scaling factor
Definition: BSKLNBaseRESPXSec2014.h:107
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::AlgFactory::GetAlgorithm
const Algorithm * GetAlgorithm(const AlgId &algid)
Definition: AlgFactory.cxx:75
genie::BSKLNBaseRESPXSec2014::fUsePauliBlocking
bool fUsePauliBlocking
account for Pauli blocking?
Definition: BSKLNBaseRESPXSec2014.h:105
genie::BSKLNBaseRESPXSec2014::fWcut
double fWcut
apply DIS/RES joining scheme &lt; Wcut
Definition: BSKLNBaseRESPXSec2014.h:99
genie::utils::res::OrbitalAngularMom
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
Definition: BaryonResUtils.cxx:531
genie::BSKLNBaseRESPXSec2014::fHAmplModelEMn
const RSHelicityAmplModelI * fHAmplModelEMn
Definition: BSKLNBaseRESPXSec2014.h:84
genie::Target::Z
int Z(void) const
Definition: Target.h:68
pINFO
#define pINFO
Definition: Messenger.h:62
genie::RSHelicityAmplModelI
Pure abstract base class. Defines the RSHelicityAmplModelI interface.
Definition: RSHelicityAmplModelI.h:27
genie::XclsTag::Resonance
Resonance_t Resonance(void) const
Definition: XclsTag.h:69
genie::Interaction::FSPrimLepton
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
Definition: Interaction.cxx:126
genie::RSHelicityAmpl::Amp20Minus
double Amp20Minus(void) const
Definition: RSHelicityAmpl.h:67
genie::ProcessInfo::IsEM
bool IsEM(void) const
Definition: ProcessInfo.cxx:193
genie::pdg::IsNeutralLepton
bool IsNeutralLepton(int pdgc)
Definition: PDGUtils.cxx:95
genie::BSKLNBaseRESPXSec2014::fGnResMaxNWidths
double fGnResMaxNWidths
limits allowed phase space for other res
Definition: BSKLNBaseRESPXSec2014.h:102
genie::FKR::C
double C
Definition: FKR.h:44
genie::utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization
double FermiMomentumForIsoscalarNucleonParametrization(const Target &target)
Definition: NuclearUtils.cxx:501
genie::BSKLNBaseRESPXSec2014::fUseRFGParametrization
bool fUseRFGParametrization
use parametrization for fermi momentum insted of table?
Definition: BSKLNBaseRESPXSec2014.h:104
genie::BSKLNBaseRESPXSec2014::fHAmplModelNCp
const RSHelicityAmplModelI * fHAmplModelNCp
Definition: BSKLNBaseRESPXSec2014.h:81
genie::FKR::Tplus
double Tplus
Definition: FKR.h:47
genie::BSKLNBaseRESPXSec2014::Integral
double Integral(const Interaction *i) const
Definition: BSKLNBaseRESPXSec2014.cxx:741
genie::AlgFactory::Instance
static AlgFactory * Instance()
Definition: AlgFactory.cxx:64
genie::FKR::B
double B
Definition: FKR.h:43
genie::FKR::Rplus
double Rplus
Definition: FKR.h:49
genie::BSKLNBaseRESPXSec2014::XSec
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
Definition: BSKLNBaseRESPXSec2014.cxx:76
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::Interaction::ExclTag
const XclsTag & ExclTag(void) const
Definition: Interaction.h:72
genie::FKR::Tminus
double Tminus
Definition: FKR.h:48
genie::pdg::IonPdgCode
int IonPdgCode(int A, int Z)
Definition: PDGUtils.cxx:71
genie::utils::kinematics::Jacobian
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
genie::RSHelicityAmpl::Amp20Plus
double Amp20Plus(void) const
Definition: RSHelicityAmpl.h:68
genie::Interaction::InitState
const InitialState & InitState(void) const
Definition: Interaction.h:69
genie::utils::res::AsString
const char * AsString(Resonance_t res)
resonance id -&gt; string
Definition: BaryonResUtils.cxx:35
genie::Interaction::ProcInfo
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
genie::BSKLNBaseRESPXSec2014::ValidProcess
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
Definition: BSKLNBaseRESPXSec2014.cxx:747
genie::FermiMomentumTable::FindClosestKF
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
Definition: FermiMomentumTable.cxx:45
genie::BSKLNBaseRESPXSec2014::fHAmplModelNCn
const RSHelicityAmplModelI * fHAmplModelNCn
Definition: BSKLNBaseRESPXSec2014.h:82
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::AlgFactory
The GENIE Algorithm Factory.
Definition: AlgFactory.h:39
genie::BSKLNBaseRESPXSec2014::fN2ResMaxNWidths
double fN2ResMaxNWidths
limits allowed phase space for n=2 res
Definition: BSKLNBaseRESPXSec2014.h:100
genie::XSecIntegratorI::Integrate
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
genie::BSKLNBaseRESPXSec2014::fZeta
double fZeta
FKR parameter Zeta.
Definition: BSKLNBaseRESPXSec2014.h:92
genie::InitialState::ProbeE
double ProbeE(RefFrame_t rf) const
Definition: InitialState.cxx:384
PDGCodes.h
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
genie::pdg::IsNegChargedLepton
bool IsNegChargedLepton(int pdgc)
Definition: PDGUtils.cxx:139
genie::FKR::S
double S
Definition: FKR.h:40
genie::BSKLNBaseRESPXSec2014::fXSecIntegrator
const XSecIntegratorI * fXSecIntegrator
Definition: BSKLNBaseRESPXSec2014.h:131
genie::BSKLNBaseRESPXSec2014::fXSecScaleEM
double fXSecScaleEM
external EM xsec scaling factor
Definition: BSKLNBaseRESPXSec2014.h:109
genie::FKR::Ta
double Ta
Definition: FKR.h:41
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
pDEBUG
#define pDEBUG
Definition: Messenger.h:63
genie::Algorithm::SubAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
Generated on Wed Jul 17 2024 14:46:38 for GENIEGenerator by   doxygen 1.8.5