GENIEGenerator
Loading...
Searching...
No Matches
genie::SmithMonizQELCCPXSec Class Reference

Computes neutrino-nucleon(nucleus) QELCC differential cross section. Is a concrete implementation of the XSecAlgorithmI interface. More...

#include <SmithMonizQELCCPXSec.h>

Inheritance diagram for genie::SmithMonizQELCCPXSec:
[legend]
Collaboration diagram for genie::SmithMonizQELCCPXSec:
[legend]

Public Member Functions

 SmithMonizQELCCPXSec ()
 SmithMonizQELCCPXSec (string config)
virtual ~SmithMonizQELCCPXSec ()
double XSec (const Interaction *i, KinePhaseSpace_t kps) const
 Compute the cross section for the input interaction.
double Integral (const Interaction *i) const
bool ValidProcess (const Interaction *i) const
 Can this cross section algorithm handle the input process?
void Configure (const Registry &config)
void Configure (string param_set)
Public Member Functions inherited from genie::XSecAlgorithmI
virtual ~XSecAlgorithmI ()
virtual bool ValidKinematics (const Interaction *i) const
 Is the input kinematical point a physically allowed one?
Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
virtual void FindConfig (void)
virtual const RegistryGetConfig (void) const
RegistryGetOwnedConfig (void)
virtual const AlgIdId (void) const
 Get algorithm ID.
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status.
virtual bool AllowReconfig (void) const
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm.
virtual void SetId (const AlgId &id)
 Set algorithm ID.
virtual void SetId (string name, string config)
const AlgorithmSubAlg (const RgKey &registry_key) const
void AdoptConfig (void)
void AdoptSubstructure (void)
virtual void Print (ostream &stream) const
 Print algorithm info.

Private Member Functions

void LoadConfig (void)
double d3sQES_dQ2dvdkF_SM (const Interaction *interaction) const
double dsQES_dQ2_SM (const Interaction *interaction) const
double d2sQES_dQ2dv_SM (const Interaction *i) const

Private Attributes

SmithMonizUtilssm_utils
double fXSecScale
 external xsec scaling factor
QELFormFactors fFormFactors
const QELFormFactorsModelIfFormFactorsModel
const XSecIntegratorIfXSecIntegrator
double fVud2
 |Vud|^2(square of magnitude ud-element of CKM-matrix)

Additional Inherited Members

Static Public Member Functions inherited from genie::Algorithm
static string BuildParamVectKey (const std::string &comm_name, unsigned int i)
static string BuildParamVectSizeKey (const std::string &comm_name)
static string BuildParamMatKey (const std::string &comm_name, unsigned int i, unsigned int j)
static string BuildParamMatRowSizeKey (const std::string &comm_name)
static string BuildParamMatColSizeKey (const std::string &comm_name)
Protected Member Functions inherited from genie::XSecAlgorithmI
 XSecAlgorithmI ()
 XSecAlgorithmI (string name)
 XSecAlgorithmI (string name, string config)
Protected Member Functions inherited from genie::Algorithm
 Algorithm ()
 Algorithm (string name)
 Algorithm (string name, string config)
void Initialize (void)
void DeleteConfig (void)
void DeleteSubstructure (void)
RegistryExtractLocalConfig (const Registry &in) const
RegistryExtractLowerConfig (const Registry &in, const string &alg_key) const
 Split an incoming configuration Registry into a block valid for the sub-algo identified by alg_key.
template<class T>
bool GetParam (const RgKey &name, T &p, bool is_top_call=true) const
template<class T>
bool GetParamDef (const RgKey &name, T &p, const T &def) const
template<class T>
int GetParamVect (const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
 Handle to load vectors of parameters.
int GetParamVectKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
template<class T>
int GetParamMat (const std::string &comm_name, TMatrixT< T > &mat, bool is_top_call=true) const
 Handle to load matrix of parameters.
template<class T>
int GetParamMatSym (const std::string &comm_name, TMatrixTSym< T > &mat, bool is_top_call=true) const
int GetParamMatKeys (const std::string &comm_name, std::vector< RgKey > &k, bool is_top_call=true) const
int AddTopRegistry (Registry *rp, bool owns=true)
 add registry with top priority, also update ownership
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership
int MergeTopRegistry (const Registry &r)
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships.
Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...)
AlgId fID
 algorithm name and configuration set
vector< Registry * > fConfVect
vector< bool > fOwnerships
 ownership for every registry in fConfVect
AlgStatus_t fStatus
 algorithm execution status
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool)

Detailed Description

Computes neutrino-nucleon(nucleus) QELCC differential cross section. Is a concrete implementation of the XSecAlgorithmI interface.

References:\n [1] R.A.Smith and E.J.Moniz,
Nuclear Physics B43, (1972) 605-622
[2] K.S. Kuzmin, V.V. Lyubushkin, V.A.Naumov, Eur. Phys. J. C54, (2008) 517-538
Author
Igor Kakorin kakor.nosp@m.in@j.nosp@m.inr.r.nosp@m.u Joint Institute for Nuclear Research
adapted from fortran code provided by:
Konstantin Kuzmin kkuzm.nosp@m.in@t.nosp@m.heor..nosp@m.jinr.nosp@m..ru Joint Institute for Nuclear Research
Vladimir Lyubushkin Joint Institute for Nuclear Research
Vadim Naumov vnaum.nosp@m.ov@t.nosp@m.heor..nosp@m.jinr.nosp@m..ru Joint Institute for Nuclear Research
based on code of:
Costas Andreopoulos <c.andreopoulos \at cern.ch> University of Liverpool
Created:\n May 05, 2017
License:\n Copyright (c) 2003-2025, The GENIE Collaboration
For the full text of the license visit http://copyright.genie-mc.org

Definition at line 52 of file SmithMonizQELCCPXSec.h.

Constructor & Destructor Documentation

◆ SmithMonizQELCCPXSec() [1/2]

SmithMonizQELCCPXSec::SmithMonizQELCCPXSec ( )

Definition at line 55 of file SmithMonizQELCCPXSec.cxx.

55 :
56XSecAlgorithmI("genie::SmithMonizQELCCPXSec")
57{
58
59}

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ SmithMonizQELCCPXSec() [2/2]

SmithMonizQELCCPXSec::SmithMonizQELCCPXSec ( string config)

Definition at line 61 of file SmithMonizQELCCPXSec.cxx.

61 :
62XSecAlgorithmI("genie::SmithMonizQELCCPXSec", config)
63{
64
65}

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ ~SmithMonizQELCCPXSec()

SmithMonizQELCCPXSec::~SmithMonizQELCCPXSec ( )
virtual

Definition at line 67 of file SmithMonizQELCCPXSec.cxx.

68{
69
70}

Member Function Documentation

◆ Configure() [1/2]

void SmithMonizQELCCPXSec::Configure ( const Registry & config)
virtual

Configure the algorithm with an external registry The registry is merged with the top level registry if it is owned, Otherwise a copy of it is added with the highest priority

Reimplemented from genie::Algorithm.

Definition at line 154 of file SmithMonizQELCCPXSec.cxx.

155{
156 Algorithm::Configure(config);
157 this->LoadConfig();
158}
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62

References genie::Algorithm::Configure(), and LoadConfig().

◆ Configure() [2/2]

void SmithMonizQELCCPXSec::Configure ( string config)
virtual

Configure the algorithm from the AlgoConfigPool based on param_set string given in input An algorithm contains a vector of registries coming from different xml configuration files, which are loaded according a very precise prioriy This methods will load a number registries in order of priority: 1) "Tunable" parameter set from CommonParametes. This is loaded with the highest prioriry and it is designed to be used for tuning procedure Usage not expected from the user. 2) For every string defined in "CommonParame" the corresponding parameter set will be loaded from CommonParameter.xml 3) parameter set specified by the config string and defined in the xml file of the algorithm 4) if config is not "Default" also the Default parameter set from the same xml file will be loaded Effectively this avoids the repetion of a parameter when it is not changed in the requested configuration

Reimplemented from genie::Algorithm.

Definition at line 160 of file SmithMonizQELCCPXSec.cxx.

161{
162 Algorithm::Configure(config);
163
164 Registry r( "SmithMonizQELCCPXSec_specific", false ) ;
165 r.Set("sm_utils_algo", RgAlg("genie::SmithMonizUtils","Default") ) ;
166
168
169 this->LoadConfig();
170}

References genie::Algorithm::Configure(), LoadConfig(), and genie::Registry::Set().

◆ d2sQES_dQ2dv_SM()

double SmithMonizQELCCPXSec::d2sQES_dQ2dv_SM ( const Interaction * i) const
private

Definition at line 311 of file SmithMonizQELCCPXSec.cxx.

312{
313 Kinematics * kinematics = interaction -> KinePtr();
314 sm_utils->SetInteraction(interaction);
315 const InitialState & init_state = interaction -> InitState();
316 // Assuming that the energy is greater of threshold.
317 // See condition in method SmithMonizQELCCXSec::Integrate
318 // interaction->InitState().ProbeE(kRfLab)<sm_utils->E_nu_thr_SM()
319 // of SmithMonizQELCCXSec.cxx
320 // if (E_nu < sm_utils->E_nu_thr_SM()) return 0;
321 // Assuming that variables Q2 and \nu are within allowable kinematic region
322 // which are specified in method: genie::utils::gsl::d2Xsec_dQ2dv::DoEval
323 double Q2 = kinematics->GetKV(kKVQ2);
324 double v = kinematics->GetKV(kKVv);
325 Range1D_t rkF = sm_utils->kFQES_SM_lim(Q2,v);
326
327 const Target & target = init_state.Tgt();
328
329
330
331// Gaussian quadratures integrate over Fermi momentum
332 double R[48]= { 0.16276744849602969579e-1,0.48812985136049731112e-1,
333 0.81297495464425558994e-1,1.13695850110665920911e-1,
334 1.45973714654896941989e-1,1.78096882367618602759e-1,
335 2.10031310460567203603e-1,2.41743156163840012328e-1,
336 2.73198812591049141487e-1,3.04364944354496353024e-1,
337 3.35208522892625422616e-1,3.65696861472313635031e-1,
338 3.95797649828908603285e-1,4.25478988407300545365e-1,
339 4.54709422167743008636e-1,4.83457973920596359768e-1,
340 5.11694177154667673586e-1,5.39388108324357436227e-1,
341 5.66510418561397168404e-1,5.93032364777572080684e-1,
342 6.18925840125468570386e-1,6.44163403784967106798e-1,
343 6.68718310043916153953e-1,6.92564536642171561344e-1,
344 7.15676812348967626225e-1,7.38030643744400132851e-1,
345 7.59602341176647498703e-1,7.80369043867433217604e-1,
346 8.00308744139140817229e-1,8.19400310737931675539e-1,
347 8.37623511228187121494e-1,8.54959033434601455463e-1,
348 8.71388505909296502874e-1,8.86894517402420416057e-1,
349 9.01460635315852341319e-1,9.15071423120898074206e-1,
350 9.27712456722308690965e-1,9.39370339752755216932e-1,
351 9.50032717784437635756e-1,9.59688291448742539300e-1,
352 9.68326828463264212174e-1,9.75939174585136466453e-1,
353 9.82517263563014677447e-1,9.88054126329623799481e-1,
354 9.92543900323762624572e-1,9.95981842987209290650e-1,
355 9.98364375863181677724e-1,9.99689503883230766828e-1};
356
357 double W[48]= { 0.00796792065552012429e-1,0.01853960788946921732e-1,
358 0.02910731817934946408e-1,0.03964554338444686674e-1,
359 0.05014202742927517693e-1,0.06058545504235961683e-1,
360 0.07096470791153865269e-1,0.08126876925698759217e-1,
361 0.09148671230783386633e-1,0.10160770535008415758e-1,
362 0.11162102099838498591e-1,0.12151604671088319635e-1,
363 0.13128229566961572637e-1,0.14090941772314860916e-1,
364 0.15038721026994938006e-1,0.15970562902562291381e-1,
365 0.16885479864245172450e-1,0.17782502316045260838e-1,
366 0.18660679627411467395e-1,0.19519081140145022410e-1,
367 0.20356797154333324595e-1,0.21172939892191298988e-1,
368 0.21966644438744349195e-1,0.22737069658329374001e-1,
369 0.23483399085926219842e-1,0.24204841792364691282e-1,
370 0.24900633222483610288e-1,0.25570036005349361499e-1,
371 0.26212340735672413913e-1,0.26826866725591762198e-1,
372 0.27412962726029242823e-1,0.27970007616848334440e-1,
373 0.28497411065085385646e-1,0.28994614150555236543e-1,
374 0.29461089958167905970e-1,0.29896344136328385984e-1,
375 0.30299915420827593794e-1,0.30671376123669149014e-1,
376 0.31010332586313837423e-1,0.31316425596861355813e-1,
377 0.31589330770727168558e-1,0.31828758894411006535e-1,
378 0.32034456231992663218e-1,0.32206204794030250669e-1,
379 0.32343822568575928429e-1,0.32447163714064269364e-1,
380 0.32516118713868835987e-1,0.32550614492363166242e-1};
381
382 double Sum = 0;
383 for(int i = 0;i<48;i++)
384 {
385 double kF = 0.5*(-R[i]*(rkF.max-rkF.min)+rkF.min+rkF.max);
386 kinematics->SetKV(kKVPn, kF);
387 Sum+=d3sQES_dQ2dvdkF_SM(interaction)*W[47-i];
388 kF = 0.5*(R[i]*(rkF.max-rkF.min)+rkF.min+rkF.max);
389 kinematics->SetKV(kKVPn, kF);
390 Sum+=d3sQES_dQ2dvdkF_SM(interaction)*W[47-i];
391 }
392
393 double xsec = 0.5*Sum*(rkF.max-rkF.min);
394
395 int nucpdgc = target.HitNucPdg();
396 int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
397
398 xsec *= NNucl; // nuclear xsec
399
400 // Apply given scaling factor
401 xsec *= fXSecScale;
402
403 return xsec;
404
405}
const Target & Tgt(void) const
double GetKV(KineVar_t kv) const
void SetKV(KineVar_t kv, double value)
double fXSecScale
external xsec scaling factor
double d3sQES_dQ2dvdkF_SM(const Interaction *interaction) const
int HitNucPdg(void) const
Definition Target.cxx:304
int N(void) const
Definition Target.h:69
int Z(void) const
Definition Target.h:68
bool IsProton(int pdgc)
Definition PDGUtils.cxx:336
double W(const Interaction *const i)
double Q2(const Interaction *const i)
@ kKVQ2
Definition KineVar.h:33
@ kKVv
Definition KineVar.h:53
@ kKVPn
Definition KineVar.h:52

References d3sQES_dQ2dvdkF_SM(), fXSecScale, genie::Target::HitNucPdg(), genie::pdg::IsProton(), genie::kKVPn, genie::kKVQ2, genie::kKVv, genie::Range1D_t::max, genie::Range1D_t::min, genie::Target::N(), sm_utils, genie::InitialState::Tgt(), and genie::Target::Z().

Referenced by XSec().

◆ d3sQES_dQ2dvdkF_SM()

double SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM ( const Interaction * interaction) const
private

Definition at line 199 of file SmithMonizQELCCPXSec.cxx.

200{
201 // Assuming that variables E_nu, Q2, \nu and kF are within allowable kinematic region
202 // which are specified in methods: genie::utils::gsl::d2Xsec_dQ2dv::DoEval and QELEventGeneratorSM::ProcessEventRecord
203 // Get kinematics & init-state parameters
204 const Kinematics & kinematics = interaction -> Kine();
205 sm_utils->SetInteraction(interaction);
206 const InitialState & init_state = interaction -> InitState();
207 const Target & target = init_state.Tgt();
208 double E_nu = init_state.ProbeE(kRfLab);
209 double Q2 = kinematics.GetKV(kKVQ2);
210 double v = kinematics.GetKV(kKVv);
211 double kF = kinematics.GetKV(kKVPn);
212 double kkF = kF*kF;
213 int nucl_pdg_ini = target.HitNucPdg();
214 int nucl_pdg_fin = genie::pdg::SwitchProtonNeutron(nucl_pdg_ini);
215
216 PDGLibrary * pdglib = PDGLibrary::Instance();
217 TParticlePDG * nucl_fin = pdglib->Find( nucl_pdg_fin );
218
219 double E_BIN = sm_utils->GetBindingEnergy();
220 double m_ini = target.HitNucMass();
221 double mm_ini = m_ini*m_ini;
222 double m_fin = nucl_fin -> Mass(); // Mass of final hadron or hadron system (GeV)
223 double mm_fin = m_fin*m_fin;
224 double m_tar = target.Mass(); // Mass of target nucleus (GeV)
225 double mm_tar = m_tar*m_tar;
226
227 // One of the xsec terms changes sign for antineutrinos
228 bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
229 int n_NT = (is_neutrino) ? +1 : -1;
230
231 double E_p = TMath::Sqrt(mm_ini+kkF)-E_BIN;
232 //|\vec{q}|
233 double qqv = v*v+Q2;
234 double qv = TMath::Sqrt(qqv);
235 double cosT_p = ((v-E_BIN)*(2*E_p+v+E_BIN)-qqv+mm_ini-mm_fin)/(2*kF*qv); //\cos\theta_p
236 if (cosT_p < -1.0 || cosT_p > 1.0 )
237 {
238 return 0.0;
239 }
240
241 double pF = TMath::Sqrt(kkF+(2*kF*qv)*cosT_p+qqv);
242
243 double E_lep = E_nu-v;
244 double m_lep = interaction->FSPrimLepton()->Mass();
245 double mm_lep = m_lep*m_lep;
246 if (E_lep < m_lep)
247 {
248 return 0.0;
249 }
250 double P_lep = TMath::Sqrt(E_lep*E_lep-mm_lep);
251 double k6 = (Q2+mm_lep)/(2*E_nu);
252 double cosT_lep= (E_lep-k6)/P_lep;
253 if (cosT_lep < -1.0 || cosT_lep > 1.0 ) return 0.0;
254
255 double cosT_k = (v+k6)/qv;
256 if (cosT_k < -1.0 || cosT_k > 1.0 ) return 0.0;
257
258 double b2_flux = (E_p-kF*cosT_k*cosT_p)*(E_p-kF*cosT_k*cosT_p);
259 double c2_flux = kkF*(1-cosT_p*cosT_p)*(1-cosT_k*cosT_k);
260
261 double k1 = fVud2*kNucleonMass2*kPi;
262 double k2 = mm_lep/(2*mm_tar);
263 double k7 = P_lep*cosT_lep;
264
265 double P_Fermi = sm_utils->GetFermiMomentum();
266 double FV_SM = 4.0*TMath::Pi()/3*TMath::Power(P_Fermi, 3);
267 double factor = k1*(m_tar*kF/(FV_SM*qv*TMath::Sqrt(b2_flux-c2_flux)))*SmithMonizUtils::rho(P_Fermi, 0.0, kF)*(1-SmithMonizUtils::rho(P_Fermi, 0.01, pF));
268
269 double a2 = kkF/kNucleonMass2;
270 double a3 = a2*cosT_p*cosT_p;
271 double a6 = kF*cosT_p/kNucleonMass;
272 double a7 = E_p/kNucleonMass;
273 double a4 = a7*a7;
274 double a5 = 2*a7*a6;
275
276 double k3 = v/qv;
277 double k4 = (3*a3-a2)/qqv;
278 double k5 = (a7-a6*k3)*m_tar/kNucleonMass;
279
280 // Calculate the QEL form factors
281 fFormFactors.Calculate(interaction);
282 double F_V = fFormFactors.F1V();
283 double F_M = fFormFactors.xiF2V();
284 double F_A = fFormFactors.FA();
285 double F_P = fFormFactors.Fp();
286 double FF_V = F_V*F_V;
287 double FF_M = F_M*F_M;
288 double FF_A = F_A*F_A;
289
290 double t = Q2/(4*kNucleonMass2);
291 double W_1 = FF_A*(1+t)+t*(F_V+F_M)*(F_V+F_M); //Ref.[1], \tilde{T}_1
292 double W_2 = FF_A+FF_V+t*FF_M; //Ref.[1], \tilde{T}_2
293 double W_3 =-2*F_A*(F_V+F_M); //Ref.[1], \tilde{T}_8
294 double W_4 =-0.5*F_V*F_M-F_A*F_P+t*F_P*F_P-0.25*(1-t)*FF_M; //Ref.[1], \tilde{T}_\alpha
295 double W_5 = FF_V+t*FF_M+FF_A;
296
297 double T_1 = 1.0*W_1+(a2-a3)*0.5*W_2; //Ref.[1], W_1
298 double T_2 = ((a2-a3)*Q2/(2*qqv)+a4-k3*(a5-k3*a3))*W_2; //Ref.[1], W_2
299 double T_3 = k5*W_3; //Ref.[1], W_8
300 double T_4 = mm_tar*(0.5*W_2*k4+1.0*W_4/kNucleonMass2+a6*W_5/(kNucleonMass*qv)); //Ref.[1], W_\alpha
301 double T_5 = k5*W_5+m_tar*(a5/qv-v*k4)*W_2;
302
303 double xsec = kGF2*factor*((E_lep-k7)*(T_1+k2*T_4)/m_tar+(E_lep+k7)*T_2/(2*m_tar)
304 +n_NT*T_3*((E_nu+E_lep)*(E_lep-k7)/(2*mm_tar)-k2)-k2*T_5)
305 *(kMw2/(kMw2+Q2))*(kMw2/(kMw2+Q2))/E_nu/kPi;
306 return xsec;
307
308
309}
int ProbePdg(void) const
double ProbeE(RefFrame_t rf) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
static PDGLibrary * Instance(void)
TParticlePDG * Find(int pdgc, bool must_exist=true)
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
static double rho(double P_Fermi, double T_Fermi, double p)
double Mass(void) const
Definition Target.cxx:224
double HitNucMass(void) const
Definition Target.cxx:233
bool IsNeutrino(int pdgc)
Definition PDGUtils.cxx:110
int SwitchProtonNeutron(int pdgc)
Definition PDGUtils.cxx:356
double Mass(Resonance_t res)
resonance mass (GeV)
@ kRfLab
Definition RefFrame.h:26

References fFormFactors, genie::PDGLibrary::Find(), genie::Interaction::FSPrimLepton(), fVud2, genie::Target::HitNucMass(), genie::Target::HitNucPdg(), genie::PDGLibrary::Instance(), genie::pdg::IsNeutrino(), genie::constants::kGF2, genie::kKVPn, genie::kKVQ2, genie::kKVv, genie::constants::kMw2, genie::constants::kNucleonMass, genie::constants::kNucleonMass2, genie::constants::kPi, genie::kRfLab, genie::Target::Mass(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::SmithMonizUtils::rho(), sm_utils, genie::pdg::SwitchProtonNeutron(), and genie::InitialState::Tgt().

Referenced by d2sQES_dQ2dv_SM(), and XSec().

◆ dsQES_dQ2_SM()

double SmithMonizQELCCPXSec::dsQES_dQ2_SM ( const Interaction * interaction) const
private

Definition at line 407 of file SmithMonizQELCCPXSec.cxx.

408{
409 // Get kinematics & init-state parameters
410 const Kinematics & kinematics = interaction -> Kine();
411 const InitialState & init_state = interaction -> InitState();
412 const Target & target = init_state.Tgt();
413
414 double E = init_state.ProbeE(kRfHitNucRest);
415 double E2 = TMath::Power(E,2);
416 double ml = interaction->FSPrimLepton()->Mass();
417 double M = target.HitNucMass();
418 double q2 = kinematics.q2();
419
420 // One of the xsec terms changes sign for antineutrinos
421 bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
422 int sign = (is_neutrino) ? -1 : 1;
423
424 // Calculate the QEL form factors
425 fFormFactors.Calculate(interaction);
426
427 double F1V = fFormFactors.F1V();
428 double xiF2V = fFormFactors.xiF2V();
429 double FA = fFormFactors.FA();
430 double Fp = fFormFactors.Fp();
431
432
433 // Calculate auxiliary parameters
434 double ml2 = TMath::Power(ml, 2);
435 double M2 = TMath::Power(M, 2);
436 double M4 = TMath::Power(M2, 2);
437 double FA2 = TMath::Power(FA, 2);
438 double Fp2 = TMath::Power(Fp, 2);
439 double F1V2 = TMath::Power(F1V, 2);
440 double xiF2V2 = TMath::Power(xiF2V, 2);
441 double Gfactor = M2*kGF2*fVud2*(kMw2/(kMw2-q2))*(kMw2/(kMw2-q2)) / (8*kPi*E2);
442 double s_u = 4*E*M + q2 - ml2;
443 double q2_M2 = q2/M2;
444
445 // Compute free nucleon differential cross section
446 double A = (0.25*(ml2-q2)/M2) * (
447 (4-q2_M2)*FA2 - (4+q2_M2)*F1V2 - q2_M2*xiF2V2*(1+0.25*q2_M2)
448 -4*q2_M2*F1V*xiF2V - (ml2/M2)*(
449 (F1V2+xiF2V2+2*F1V*xiF2V)+(FA2+4*Fp2+4*FA*Fp)+(q2_M2-4)*Fp2));
450 double B = -1 * q2_M2 * FA*(F1V+xiF2V);
451 double C = 0.25*(FA2 + F1V2 - 0.25*q2_M2*xiF2V2);
452
453 double xsec = Gfactor * (A + sign*B*s_u/M2 + C*s_u*s_u/M4);
454
455 // Apply given scaling factor
456 xsec *= fXSecScale;
457
458 // Pauli-correction factor for deuterium, we formally apply this factor for He-3 and tritium,
459 // because RFG model is not applicable for them.
460 if (1<target.A() && target.A()<4)
461 {
462 double Q2 = -q2;
463 double fQES_Pauli = 1.0-0.529*TMath::Exp((Q2*(228.0-531.0*Q2)-48.0)*Q2);
464 xsec *= fQES_Pauli;
465 }
466
467 int nucpdgc = target.HitNucPdg();
468 int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
469
470 xsec *= NNucl; // nuclear xsec
471
472 // Apply radiative correction to the cross section for IBD processes
473 // Refs:
474 // 1) I.S. Towner, Phys. Rev. C 58 (1998) 1288;
475 // 2) J.F. Beacom, S.J. Parke, Phys. Rev. D 64 (2001) 091302;
476 // 3) A. Kurylov, M.J. Ramsey-Musolf, P. Vogel, Phys. Rev. C 65 (2002) 055501;
477 // 4) A. Kurylov, M.J. Ramsey-Musolf, P. Vogel, Phys. Rev. C 67 (2003) 035502.
478 double rc = 1.0;
479 if ( (target.IsProton() && pdg::IsAntiNuE(init_state.ProbePdg())) || (target.IsNeutron() && pdg::IsNuE(init_state.ProbePdg()) ))
480 {
481 const double mp = kProtonMass;
482 const double mp2 = kProtonMass2;
483 const double mn2 = kNeutronMass2;
484 const double Ee = E + ( (q2 - mn2 + mp2) / 2.0 / mp );
485 assert(Ee > 0.0); // must be non-zero and positive
486 rc = 6.0 + (1.5 * TMath::Log(kProtonMass / 2.0 / Ee));
487 rc += 1.2 * TMath::Power((kElectronMass / Ee), 1.5);
488 rc *= kAem / kPi;
489 rc += 1.0;
490 }
491
492 xsec *= rc;
493 return xsec;
494}
double q2(bool selected=false) const
bool IsNeutron(void) const
Definition Target.cxx:267
int A(void) const
Definition Target.h:70
bool IsProton(void) const
Definition Target.cxx:262
bool IsNuE(int pdgc)
Definition PDGUtils.cxx:158
bool IsAntiNuE(int pdgc)
Definition PDGUtils.cxx:173
@ kRfHitNucRest
Definition RefFrame.h:30

References genie::Target::A(), fFormFactors, genie::Interaction::FSPrimLepton(), fVud2, fXSecScale, genie::Target::HitNucMass(), genie::Target::HitNucPdg(), genie::pdg::IsAntiNuE(), genie::pdg::IsNeutrino(), genie::Target::IsNeutron(), genie::pdg::IsNuE(), genie::pdg::IsProton(), genie::Target::IsProton(), genie::constants::kAem, genie::constants::kElectronMass, genie::constants::kGF2, genie::constants::kMw2, genie::constants::kNeutronMass2, genie::constants::kPi, genie::constants::kProtonMass, genie::constants::kProtonMass2, genie::kRfHitNucRest, genie::Target::N(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::InitialState::Tgt(), and genie::Target::Z().

Referenced by XSec().

◆ Integral()

double SmithMonizQELCCPXSec::Integral ( const Interaction * i) const
virtual

Integrate the model over the kinematic phase space available to the input interaction (kinematical cuts can be included)

Implements genie::XSecAlgorithmI.

Definition at line 125 of file SmithMonizQELCCPXSec.cxx.

126{
127 return fXSecIntegrator->Integrate(this,in);
128
129}
const XSecIntegratorI * fXSecIntegrator

References fXSecIntegrator.

◆ LoadConfig()

void SmithMonizQELCCPXSec::LoadConfig ( void )
private

Definition at line 172 of file SmithMonizQELCCPXSec.cxx.

173{
174
175 // Cross section scaling factor
176 GetParamDef( "QEL-CC-XSecScale", fXSecScale, 1. ) ;
177
178 double Vud;
179 GetParam( "CKM-Vud", Vud ) ;
180 fVud2 = TMath::Power( Vud, 2 );
181
182 // load QEL form factors model
183 fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
184 this->SubAlg("FormFactorsAlg"));
185 assert(fFormFactorsModel);
186 fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
187
188 // load XSec Integrators
190 dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
191 assert(fXSecIntegrator);
192
193 sm_utils = const_cast<genie::SmithMonizUtils *>(
194 dynamic_cast<const genie::SmithMonizUtils *>(
195 this -> SubAlg( "sm_utils_algo" ) ) ) ;
196
197}
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
bool GetParamDef(const RgKey &name, T &p, const T &def) const
const Algorithm * SubAlg(const RgKey &registry_key) const
const QELFormFactorsModelI * fFormFactorsModel

References fFormFactors, fFormFactorsModel, fVud2, fXSecIntegrator, fXSecScale, genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), sm_utils, and genie::Algorithm::SubAlg().

Referenced by Configure(), and Configure().

◆ ValidProcess()

bool SmithMonizQELCCPXSec::ValidProcess ( const Interaction * i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 131 of file SmithMonizQELCCPXSec.cxx.

132{
133 if(interaction->TestBit(kISkipProcessChk)) return true;
134
135 const InitialState & init_state = interaction->InitState();
136 const ProcessInfo & proc_info = interaction->ProcInfo();
137
138 if(!proc_info.IsQuasiElastic()) return false;
139
140 int nuc = init_state.Tgt().HitNucPdg();
141 int nu = init_state.ProbePdg();
142
143 bool isP = pdg::IsProton(nuc);
144 bool isN = pdg::IsNeutron(nuc);
145 bool isnu = pdg::IsNeutrino(nu);
146 bool isnub = pdg::IsAntiNeutrino(nu);
147
148 bool prcok = proc_info.IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
149 if(!prcok) return false;
150
151 return true;
152}
bool IsWeakCC(void) const
bool IsQuasiElastic(void) const
bool IsNeutron(int pdgc)
Definition PDGUtils.cxx:341
bool IsAntiNeutrino(int pdgc)
Definition PDGUtils.cxx:118
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition Interaction.h:47

References genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IsAntiNeutrino(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::ProcessInfo::IsQuasiElastic(), genie::ProcessInfo::IsWeakCC(), genie::kISkipProcessChk, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by XSec().

◆ XSec()

double SmithMonizQELCCPXSec::XSec ( const Interaction * i,
KinePhaseSpace_t k ) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 72 of file SmithMonizQELCCPXSec.cxx.

74{
75 double xsec = 0. ;
76 // dimension of kine phase space
77 std::string s = KinePhaseSpace::AsString(kps);
78 int kpsdim = s!="<|E>"?1 + std::count(s.begin(), s.begin()+s.find('}'), ','):0;
79
80 if(!this -> ValidProcess (interaction) )
81 {
82 LOG("SmithMoniz",pWARN) << "not a valid process";
83 return 0.;
84 }
85
86 if(kpsdim == 1)
87 {
88 if(! this -> ValidKinematics (interaction) )
89 {
90 LOG("SmithMoniz",pWARN) << "not valid kinematics";
91 return 0.;
92 }
93 xsec = this->dsQES_dQ2_SM(interaction);
94 }
95
96 if(kpsdim == 2)
97 {
98 xsec = this->d2sQES_dQ2dv_SM(interaction);
99 }
100
101 if(kpsdim == 3)
102 {
103 xsec = this->d3sQES_dQ2dvdkF_SM(interaction);
104 }
105
106
107 // The algorithm computes d^1xsec/dQ2, d^2xsec/dQ2dv or d^3xsec/dQ2dvdp
108 // Check whether variable tranformation is needed
109 if ( kps != kPSQ2fE && kps != kPSQ2vfE )
110 {
111 double J = 1.;
112 if (kpsdim == 1)
113 J = utils::kinematics::Jacobian(interaction, kPSQ2fE, kps);
114 else if (kpsdim == 2)
115 J = utils::kinematics::Jacobian(interaction, kPSQ2vfE, kps);
116 else if (kpsdim == 3)
117 J = utils::kinematics::Jacobian(interaction, kPSQ2vpfE, kps);
118 xsec *= J;
119 }
120
121 return xsec;
122
123}
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
#define pWARN
Definition Messenger.h:60
static string AsString(KinePhaseSpace_t kps)
double d2sQES_dQ2dv_SM(const Interaction *i) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double dsQES_dQ2_SM(const Interaction *interaction) const
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
double J(double q0, double q3, double Enu, double ml)
Definition MECUtils.cxx:147

References genie::KinePhaseSpace::AsString(), d2sQES_dQ2dv_SM(), d3sQES_dQ2dvdkF_SM(), dsQES_dQ2_SM(), genie::utils::kinematics::Jacobian(), genie::kPSQ2fE, genie::kPSQ2vfE, genie::kPSQ2vpfE, LOG, pWARN, genie::XSecAlgorithmI::ValidKinematics(), and ValidProcess().

Member Data Documentation

◆ fFormFactors

QELFormFactors genie::SmithMonizQELCCPXSec::fFormFactors
mutableprivate

Definition at line 78 of file SmithMonizQELCCPXSec.h.

Referenced by d3sQES_dQ2dvdkF_SM(), dsQES_dQ2_SM(), and LoadConfig().

◆ fFormFactorsModel

const QELFormFactorsModelI* genie::SmithMonizQELCCPXSec::fFormFactorsModel
private

Definition at line 79 of file SmithMonizQELCCPXSec.h.

Referenced by LoadConfig().

◆ fVud2

double genie::SmithMonizQELCCPXSec::fVud2
private

|Vud|^2(square of magnitude ud-element of CKM-matrix)

Definition at line 81 of file SmithMonizQELCCPXSec.h.

Referenced by d3sQES_dQ2dvdkF_SM(), dsQES_dQ2_SM(), and LoadConfig().

◆ fXSecIntegrator

const XSecIntegratorI* genie::SmithMonizQELCCPXSec::fXSecIntegrator
private

Definition at line 80 of file SmithMonizQELCCPXSec.h.

Referenced by Integral(), and LoadConfig().

◆ fXSecScale

double genie::SmithMonizQELCCPXSec::fXSecScale
private

external xsec scaling factor

Definition at line 77 of file SmithMonizQELCCPXSec.h.

Referenced by d2sQES_dQ2dv_SM(), dsQES_dQ2_SM(), and LoadConfig().

◆ sm_utils

SmithMonizUtils* genie::SmithMonizQELCCPXSec::sm_utils
mutableprivate

Definition at line 70 of file SmithMonizQELCCPXSec.h.

Referenced by d2sQES_dQ2dv_SM(), d3sQES_dQ2dvdkF_SM(), and LoadConfig().


The documentation for this class was generated from the following files: