Public Member Functions | Private Member Functions | Private Attributes | List of all members
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:
genie::XSecAlgorithmI genie::Algorithm

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. More...
 
double Integral (const Interaction *i) const
 
bool ValidProcess (const Interaction *i) const
 Can this cross section algorithm handle the input process? More...
 
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? More...
 
- 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. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
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. More...
 

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 More...
 
QELFormFactors fFormFactors
 
const QELFormFactorsModelIfFormFactorsModel
 
const XSecIntegratorIfXSecIntegrator
 
double fVud2
 |Vud|^2(square of magnitude ud-element of CKM-matrix) More...
 
int fn_NT
 
double fQ2
 
double fv
 
double fE_nu
 
double fE_lep
 
double fmm_ini
 
double fmm_fin
 
double fm_tar
 
double fmm_tar
 
double fk1
 
double fk2
 
double fk7
 
double fqv
 
double fqqv
 
double fcosT_k
 
double fF_V
 
double fF_M
 
double fF_A
 
double fF_P
 
double fFF_V
 
double fFF_M
 
double fFF_A
 
double fW_1
 
double fW_2
 
double fW_3
 
double fW_4
 
double fW_5
 

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)
 
- 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. More...
 
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. More...
 
int GetParamVectKeys (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 More...
 
int AddLowRegistry (Registry *rp, bool owns=true)
 add registry with lowest priority, also update ownership More...
 
int MergeTopRegistry (const Registry &r)
 
int AddTopRegisties (const vector< Registry * > &rs, bool owns=false)
 Add registries with top priority, also udated Ownerships. More...
 
- Protected Attributes inherited from genie::Algorithm
bool fAllowReconfig
 
bool fOwnsSubstruc
 true if it owns its substructure (sub-algs,...) More...
 
AlgId fID
 algorithm name and configuration set More...
 
vector< Registry * > fConfVect
 
vector< boolfOwnerships
 ownership for every registry in fConfVect More...
 
AlgStatus_t fStatus
 algorithm execution status More...
 
AlgMapfOwnedSubAlgMp
 local pool for owned sub-algs (taken out of the factory pool) More...
 

Detailed Description

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

[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 <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

May 05, 2017

Copyright (c) 2003-2020, 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::SmithMonizQELCCPXSec ( )

Definition at line 55 of file SmithMonizQELCCPXSec.cxx.

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

Definition at line 61 of file SmithMonizQELCCPXSec.cxx.

61  :
62 XSecAlgorithmI("genie::SmithMonizQELCCPXSec", config)
63 {
64 
65 }
static Config * config
Definition: config.cpp:1054
SmithMonizQELCCPXSec::~SmithMonizQELCCPXSec ( )
virtual

Definition at line 67 of file SmithMonizQELCCPXSec.cxx.

68 {
69 
70 }

Member Function Documentation

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 148 of file SmithMonizQELCCPXSec.cxx.

149 {
150  Algorithm::Configure(config);
151  this->LoadConfig();
152 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
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 154 of file SmithMonizQELCCPXSec.cxx.

155 {
157 
158  Registry r( "SmithMonizQELCCPXSec_specific", false ) ;
159  r.Set("sm_utils_algo", RgAlg("genie::SmithMonizUtils","Default") ) ;
160 
162 
163  this->LoadConfig();
164 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
double SmithMonizQELCCPXSec::d2sQES_dQ2dv_SM ( const Interaction i) const
private

Definition at line 240 of file SmithMonizQELCCPXSec.cxx.

241 {
242  Kinematics * kinematics = interaction -> KinePtr();
244  const InitialState & init_state = interaction -> InitState();
245  fE_nu = init_state.ProbeE(kRfLab);
246  if (fE_nu < sm_utils->E_nu_thr_SM()) return 0;
247  fQ2 = kinematics->GetKV(kKVQ2);
248  fv = kinematics->GetKV(kKVv);
250 
251  const Target & target = init_state.Tgt();
252  PDGLibrary * pdglib = PDGLibrary::Instance();
253 
254  // One of the xsec terms changes sign for antineutrinos
255  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
256  fn_NT = (is_neutrino) ? +1 : -1;
257 
258  int nucl_pdg_ini = target.HitNucPdg();
259  double m_ini = target.HitNucMass();
260  fmm_ini = TMath::Power(m_ini, 2);
261  int nucl_pdg_fin = genie::pdg::SwitchProtonNeutron(nucl_pdg_ini);
262  TParticlePDG * nucl_fin = pdglib->Find( nucl_pdg_fin );
263  double m_fin = nucl_fin -> Mass(); // Mass of final hadron or hadron system (GeV)
264  fmm_fin = TMath::Power(m_fin, 2);
265  fm_tar = target.Mass(); // Mass of target nucleus (GeV)
266  fmm_tar = TMath::Power(fm_tar, 2);
267 
268  fE_lep = fE_nu-fv;
269  double m_lep = interaction->FSPrimLepton()->Mass();
270  double mm_lep = m_lep*m_lep;
271  if (fE_lep < m_lep) return 0.0;
272  double P_lep = TMath::Sqrt(fE_lep*fE_lep-mm_lep);
273  double k6 = (fQ2+mm_lep)/(2*fE_nu);
274  double cosT_lep= (fE_lep-k6)/P_lep;
275  if (cosT_lep < -1.0 || cosT_lep > 1.0 ) return 0.0;
276  //|\vec{q}|
277  fqqv = fv*fv+fQ2;
278  fqv = TMath::Sqrt(fqqv);
279  fcosT_k = (fv+k6)/fqv;
280  if (fcosT_k < -1.0 || fcosT_k > 1.0 ) return 0.0;
281 
283  fk2 = mm_lep/(2*fmm_tar);
284  fk7 = P_lep*cosT_lep;
285 
286 
287  // Calculate the QEL form factors
289  fF_V = fFormFactors.F1V();
290  fF_M = fFormFactors.xiF2V();
291  fF_A = fFormFactors.FA();
292  fF_P = fFormFactors.Fp();
293  fFF_V = fF_V*fF_V;
294  fFF_M = fF_M*fF_M;
295  fFF_A = fF_A*fF_A;
296 
297  double t = fQ2/(4*kNucleonMass2);
298  fW_1 = fFF_A*(1+t)+t*(fF_V+fF_M)*(fF_V+fF_M); //Ref.[1], \tilde{T}_1
299  fW_2 = fFF_A+fFF_V+t*fFF_M; //Ref.[1], \tilde{T}_2
300  fW_3 =-2*fF_A*(fF_V+fF_M); //Ref.[1], \tilde{T}_8
301  fW_4 =-0.5*fF_V*fF_M-fF_A*fF_P+t*fF_P*fF_P-0.25*(1-t)*fFF_M; //Ref.[1], \tilde{T}_\alpha
302  fW_5 = fFF_V+t*fFF_M+fFF_A;
303 
304 // Gaussian quadratures integrate over Fermi momentum
305  double R[48]= { 0.16276744849602969579e-1,0.48812985136049731112e-1,
306  0.81297495464425558994e-1,1.13695850110665920911e-1,
307  1.45973714654896941989e-1,1.78096882367618602759e-1,
308  2.10031310460567203603e-1,2.41743156163840012328e-1,
309  2.73198812591049141487e-1,3.04364944354496353024e-1,
310  3.35208522892625422616e-1,3.65696861472313635031e-1,
311  3.95797649828908603285e-1,4.25478988407300545365e-1,
312  4.54709422167743008636e-1,4.83457973920596359768e-1,
313  5.11694177154667673586e-1,5.39388108324357436227e-1,
314  5.66510418561397168404e-1,5.93032364777572080684e-1,
315  6.18925840125468570386e-1,6.44163403784967106798e-1,
316  6.68718310043916153953e-1,6.92564536642171561344e-1,
317  7.15676812348967626225e-1,7.38030643744400132851e-1,
318  7.59602341176647498703e-1,7.80369043867433217604e-1,
319  8.00308744139140817229e-1,8.19400310737931675539e-1,
320  8.37623511228187121494e-1,8.54959033434601455463e-1,
321  8.71388505909296502874e-1,8.86894517402420416057e-1,
322  9.01460635315852341319e-1,9.15071423120898074206e-1,
323  9.27712456722308690965e-1,9.39370339752755216932e-1,
324  9.50032717784437635756e-1,9.59688291448742539300e-1,
325  9.68326828463264212174e-1,9.75939174585136466453e-1,
326  9.82517263563014677447e-1,9.88054126329623799481e-1,
327  9.92543900323762624572e-1,9.95981842987209290650e-1,
328  9.98364375863181677724e-1,9.99689503883230766828e-1};
329 
330  double W[48]= { 0.00796792065552012429e-1,0.01853960788946921732e-1,
331  0.02910731817934946408e-1,0.03964554338444686674e-1,
332  0.05014202742927517693e-1,0.06058545504235961683e-1,
333  0.07096470791153865269e-1,0.08126876925698759217e-1,
334  0.09148671230783386633e-1,0.10160770535008415758e-1,
335  0.11162102099838498591e-1,0.12151604671088319635e-1,
336  0.13128229566961572637e-1,0.14090941772314860916e-1,
337  0.15038721026994938006e-1,0.15970562902562291381e-1,
338  0.16885479864245172450e-1,0.17782502316045260838e-1,
339  0.18660679627411467395e-1,0.19519081140145022410e-1,
340  0.20356797154333324595e-1,0.21172939892191298988e-1,
341  0.21966644438744349195e-1,0.22737069658329374001e-1,
342  0.23483399085926219842e-1,0.24204841792364691282e-1,
343  0.24900633222483610288e-1,0.25570036005349361499e-1,
344  0.26212340735672413913e-1,0.26826866725591762198e-1,
345  0.27412962726029242823e-1,0.27970007616848334440e-1,
346  0.28497411065085385646e-1,0.28994614150555236543e-1,
347  0.29461089958167905970e-1,0.29896344136328385984e-1,
348  0.30299915420827593794e-1,0.30671376123669149014e-1,
349  0.31010332586313837423e-1,0.31316425596861355813e-1,
350  0.31589330770727168558e-1,0.31828758894411006535e-1,
351  0.32034456231992663218e-1,0.32206204794030250669e-1,
352  0.32343822568575928429e-1,0.32447163714064269364e-1,
353  0.32516118713868835987e-1,0.32550614492363166242e-1};
354 
355  double Sum = 0;
356  for(int i = 0;i<48;i++)
357  {
358  double kF = 0.5*(-R[i]*(rkF.max-rkF.min)+rkF.min+rkF.max);
359  kinematics->SetKV(kKVPn, kF);
360  Sum+=d3sQES_dQ2dvdkF_SM(interaction)*W[47-i];
361  kF = 0.5*(R[i]*(rkF.max-rkF.min)+rkF.min+rkF.max);
362  kinematics->SetKV(kKVPn, kF);
363  Sum+=d3sQES_dQ2dvdkF_SM(interaction)*W[47-i];
364  }
365 
366  double xsec = 0.5*Sum*(rkF.max-rkF.min);
367 
368  int nucpdgc = target.HitNucPdg();
369  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
370 
371  xsec *= NNucl; // nuclear xsec
372 
373  // Apply given scaling factor
374  xsec *= fXSecScale;
375 
376  return xsec;
377 
378 }
void SetInteraction(const Interaction *i)
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
int HitNucPdg(void) const
Definition: Target.cxx:304
double fXSecScale
external xsec scaling factor
A simple [min,max] interval for doubles.
Definition: Range1.h:42
double HitNucMass(void) const
Definition: Target.cxx:233
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
int SwitchProtonNeutron(int pdgc)
Definition: PDGUtils.cxx:353
double d3sQES_dQ2dvdkF_SM(const Interaction *interaction) const
double Mass(Resonance_t res)
resonance mass (GeV)
double Mass(void) const
Definition: Target.cxx:224
Range1D_t kFQES_SM_lim(double nu, double Q2) const
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
int Z(void) const
Definition: Target.h:68
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
static const double kNucleonMass2
Definition: Constants.h:89
double xiF2V(void) const
Get the computed form factor xi*F2V.
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
void SetKV(KineVar_t kv, double value)
Definition: Kinematics.cxx:335
double max
Definition: Range1.h:53
int N(void) const
Definition: Target.h:69
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:57
Singleton class to load & serve a TDatabasePDG.
Definition: PDGLibrary.h:32
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
double min
Definition: Range1.h:52
double F1V(void) const
Get the computed form factor F1V.
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:75
const Target & Tgt(void) const
Definition: InitialState.h:66
double Fp(void) const
Get the computed form factor Fp.
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:37
double FA(void) const
Get the computed form factor FA.
Initial State information.
Definition: InitialState.h:48
double SmithMonizQELCCPXSec::d3sQES_dQ2dvdkF_SM ( const Interaction interaction) const
private

Definition at line 193 of file SmithMonizQELCCPXSec.cxx.

194 {
195  // Get kinematics & init-state parameters
196  const Kinematics & kinematics = interaction -> Kine();
197 
198  double kF = kinematics.GetKV(kKVPn);
199  double kkF = kF*kF;
200  double P_Fermi, E_nuBIN;
201 
202  E_nuBIN = sm_utils->GetBindingEnergy();
203 
204  double E_p = TMath::Sqrt(fmm_ini+kkF)-E_nuBIN;
205  double cosT_p = ((fv-E_nuBIN)*(2*E_p+fv+E_nuBIN)-fqqv+fmm_ini-fmm_fin)/(2*kF*fqv); //\cos\theta_p
206  if (cosT_p < -1.0 || cosT_p > 1.0 ) return 0.0;
207  double pF = TMath::Sqrt(kkF+(2*kF*fqv)*cosT_p+fqqv);
208  double b2_flux = (E_p-kF*fcosT_k*cosT_p)*(E_p-kF*fcosT_k*cosT_p);
209  double c2_flux = kkF*(1-cosT_p*cosT_p)*(1-fcosT_k*fcosT_k);
210 
211  P_Fermi = sm_utils->GetFermiMomentum();
212  double FV_SM = 4.0*TMath::Pi()/3*TMath::Power(P_Fermi, 3);
213  double factor = fk1*(fm_tar*kF/(FV_SM*fqv*TMath::Sqrt(b2_flux-c2_flux)))*SmithMonizUtils::rho(P_Fermi, 0.0, kF)*(1-SmithMonizUtils::rho(P_Fermi, 0.01, pF));
214 
215  double a2 = kkF/kNucleonMass2;
216  double a3 = a2*cosT_p*cosT_p;
217  double a6 = kF*cosT_p/kNucleonMass;
218  double a7 = E_p/kNucleonMass;
219  double a4 = a7*a7;
220  double a5 = 2*a7*a6;
221 
222  double k3 = fv/fqv;
223  double k4 = (3*a3-a2)/fqqv;
224  double k5 = (a7-a6*k3)*fm_tar/kNucleonMass;
225 
226  double T_1 = 1.0*fW_1+(a2-a3)*0.5*fW_2; //Ref.[1], W_1
227  double T_2 = ((a2-a3)*fQ2/(2*fqqv)+a4-k3*(a5-k3*a3))*fW_2; //Ref.[1], W_2
228  double T_3 = k5*fW_3; //Ref.[1], W_8
229  double T_4 = fmm_tar*(0.5*fW_2*k4+1.0*fW_4/kNucleonMass2+a6*fW_5/(kNucleonMass*fqv)); //Ref.[1], W_\alpha
230  double T_5 = k5*fW_5+fm_tar*(a5/fqv-fv*k4)*fW_2;
231 
232  double xsec = kGF2*factor*((fE_lep-fk7)*(T_1+fk2*T_4)/fm_tar+(fE_lep+fk7)*T_2/(2*fm_tar)
233  +fn_NT*T_3*((fE_nu+fE_lep)*(fE_lep-fk7)/(2*fmm_tar)-fk2)-fk2*T_5)
234  *(kMw2/(kMw2+fQ2))*(kMw2/(kMw2+fQ2))/fE_nu/kPi;
235  return xsec;
236 
237 
238 }
#define a6
static const double kMw2
Definition: Constants.h:93
static const double kNucleonMass
Definition: Constants.h:77
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double GetBindingEnergy(void) const
#define a5
#define a2
#define a3
double GetKV(KineVar_t kv) const
Definition: Kinematics.cxx:323
static const double kNucleonMass2
Definition: Constants.h:89
static double rho(double P_Fermi, double T_Fermi, double p)
#define a4
static const double kGF2
Definition: Constants.h:59
static const double kPi
Definition: Constants.h:37
double GetFermiMomentum(void) const
double SmithMonizQELCCPXSec::dsQES_dQ2_SM ( const Interaction interaction) const
private

Definition at line 380 of file SmithMonizQELCCPXSec.cxx.

381 {
382  // Get kinematics & init-state parameters
383  const Kinematics & kinematics = interaction -> Kine();
384  const InitialState & init_state = interaction -> InitState();
385  const Target & target = init_state.Tgt();
386 
387  double E = init_state.ProbeE(kRfHitNucRest);
388  double E2 = TMath::Power(E,2);
389  double ml = interaction->FSPrimLepton()->Mass();
390  double M = target.HitNucMass();
391  double q2 = kinematics.q2();
392 
393  // One of the xsec terms changes sign for antineutrinos
394  bool is_neutrino = pdg::IsNeutrino(init_state.ProbePdg());
395  int sign = (is_neutrino) ? -1 : 1;
396 
397  // Calculate the QEL form factors
398  fFormFactors.Calculate(interaction);
399 
400  double F1V = fFormFactors.F1V();
401  double xiF2V = fFormFactors.xiF2V();
402  double FA = fFormFactors.FA();
403  double Fp = fFormFactors.Fp();
404 
405 
406  // Calculate auxiliary parameters
407  double ml2 = TMath::Power(ml, 2);
408  double M2 = TMath::Power(M, 2);
409  double M4 = TMath::Power(M2, 2);
410  double FA2 = TMath::Power(FA, 2);
411  double Fp2 = TMath::Power(Fp, 2);
412  double F1V2 = TMath::Power(F1V, 2);
413  double xiF2V2 = TMath::Power(xiF2V, 2);
414  double Gfactor = M2*kGF2*fVud2*(kMw2/(kMw2-q2))*(kMw2/(kMw2-q2)) / (8*kPi*E2);
415  double s_u = 4*E*M + q2 - ml2;
416  double q2_M2 = q2/M2;
417 
418  // Compute free nucleon differential cross section
419  double A = (0.25*(ml2-q2)/M2) * (
420  (4-q2_M2)*FA2 - (4+q2_M2)*F1V2 - q2_M2*xiF2V2*(1+0.25*q2_M2)
421  -4*q2_M2*F1V*xiF2V - (ml2/M2)*(
422  (F1V2+xiF2V2+2*F1V*xiF2V)+(FA2+4*Fp2+4*FA*Fp)+(q2_M2-4)*Fp2));
423  double B = -1 * q2_M2 * FA*(F1V+xiF2V);
424  double C = 0.25*(FA2 + F1V2 - 0.25*q2_M2*xiF2V2);
425 
426  double xsec = Gfactor * (A + sign*B*s_u/M2 + C*s_u*s_u/M4);
427 
428  // Apply given scaling factor
429  xsec *= fXSecScale;
430 
431  // Deuterium and tritium is a special case
432  if (target.A()>1 && target.A()<4)
433  {
434  double Q2 = -q2;
435  double fQES_Pauli = 1.0-0.529*TMath::Exp((Q2*(228.0-531.0*Q2)-48.0)*Q2);
436  xsec *= fQES_Pauli;
437  }
438 
439  int nucpdgc = target.HitNucPdg();
440  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
441 
442  xsec *= NNucl; // nuclear xsec
443 
444  // Apply radiative correction to the cross section for IBD processes
445  // Refs:
446  // 1) I.S. Towner, Phys. Rev. C 58 (1998) 1288;
447  // 2) J.F. Beacom, S.J. Parke, Phys. Rev. D 64 (2001) 091302;
448  // 3) A. Kurylov, M.J. Ramsey-Musolf, P. Vogel, Phys. Rev. C 65 (2002) 055501;
449  // 4) A. Kurylov, M.J. Ramsey-Musolf, P. Vogel, Phys. Rev. C 67 (2003) 035502.
450  double rc = 1.0;
451  if ( (target.IsProton() && pdg::IsAntiNuE(init_state.ProbePdg())) || (target.IsNeutron() && pdg::IsNuE(init_state.ProbePdg()) ))
452  {
453  const double mp = kProtonMass;
454  const double mp2 = kProtonMass2;
455  const double mn2 = kNeutronMass2;
456  const double Ee = E + ( (q2 - mn2 + mp2) / 2.0 / mp );
457  assert(Ee > 0.0); // must be non-zero and positive
458  rc = 6.0 + (1.5 * TMath::Log(kProtonMass / 2.0 / Ee));
459  rc += 1.2 * TMath::Power((kElectronMass / Ee), 1.5);
460  rc *= kAem / kPi;
461  rc += 1.0;
462  }
463 
464  xsec *= rc;
465  return xsec;
466 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
static const double kMw2
Definition: Constants.h:93
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
int HitNucPdg(void) const
Definition: Target.cxx:304
bool IsNeutron(void) const
Definition: Target.cxx:267
double fXSecScale
external xsec scaling factor
int A(void) const
Definition: Target.h:70
double HitNucMass(void) const
Definition: Target.cxx:233
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
bool IsNuE(int pdgc)
Definition: PDGUtils.cxx:155
static const double kElectronMass
Definition: Constants.h:70
static const double kAem
Definition: Constants.h:56
double q2(bool selected=false) const
Definition: Kinematics.cxx:141
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int ProbePdg(void) const
Definition: InitialState.h:64
int Z(void) const
Definition: Target.h:68
double xiF2V(void) const
Get the computed form factor xi*F2V.
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
void Calculate(const Interaction *interaction)
Compute the form factors for the input interaction using the attached model.
int sign(double val)
Definition: UtilFunc.cxx:104
int N(void) const
Definition: Target.h:69
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
Definition: 018_def.c:13
#define A
Definition: memgrp.cpp:38
double F1V(void) const
Get the computed form factor F1V.
const Target & Tgt(void) const
Definition: InitialState.h:66
static const double kProtonMass2
Definition: Constants.h:87
static const double kGF2
Definition: Constants.h:59
double Fp(void) const
Get the computed form factor Fp.
static const double kNeutronMass2
Definition: Constants.h:88
bool IsProton(void) const
Definition: Target.cxx:262
double ProbeE(RefFrame_t rf) const
static const double kPi
Definition: Constants.h:37
double FA(void) const
Get the computed form factor FA.
bool IsAntiNuE(int pdgc)
Definition: PDGUtils.cxx:170
static const double kProtonMass
Definition: Constants.h:75
Initial State information.
Definition: InitialState.h:48
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 119 of file SmithMonizQELCCPXSec.cxx.

120 {
121  return fXSecIntegrator->Integrate(this,in);
122 
123 }
const XSecIntegratorI * fXSecIntegrator
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
void SmithMonizQELCCPXSec::LoadConfig ( void  )
private

Definition at line 166 of file SmithMonizQELCCPXSec.cxx.

167 {
168 
169  // Cross section scaling factor
170  GetParamDef( "QEL-CC-XSecScale", fXSecScale, 1. ) ;
171 
172  double Vud;
173  GetParam( "CKM-Vud", Vud ) ;
174  fVud2 = TMath::Power( Vud, 2 );
175 
176  // load QEL form factors model
177  fFormFactorsModel = dynamic_cast<const QELFormFactorsModelI *> (
178  this->SubAlg("FormFactorsAlg"));
179  assert(fFormFactorsModel);
180  fFormFactors.SetModel(fFormFactorsModel); // <-- attach algorithm
181 
182  // load XSec Integrators
184  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator"));
185  assert(fXSecIntegrator);
186 
187  sm_utils = const_cast<genie::SmithMonizUtils *>(
188  dynamic_cast<const genie::SmithMonizUtils *>(
189  this -> SubAlg( "sm_utils_algo" ) ) ) ;
190 
191 }
Cross Section Integrator Interface.
double fXSecScale
external xsec scaling factor
void SetModel(const QELFormFactorsModelI *model)
Attach an algorithm.
Pure abstract base class. Defines the QELFormFactorsModelI interface to be implemented by any algorit...
const QELFormFactorsModelI * fFormFactorsModel
const XSecIntegratorI * fXSecIntegrator
Contains auxiliary functions for Smith-Moniz model. .
double fVud2
|Vud|^2(square of magnitude ud-element of CKM-matrix)
bool GetParamDef(const RgKey &name, T &p, const T &def) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
bool SmithMonizQELCCPXSec::ValidProcess ( const Interaction i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 125 of file SmithMonizQELCCPXSec.cxx.

126 {
127  if(interaction->TestBit(kISkipProcessChk)) return true;
128 
129  const InitialState & init_state = interaction->InitState();
130  const ProcessInfo & proc_info = interaction->ProcInfo();
131 
132  if(!proc_info.IsQuasiElastic()) return false;
133 
134  int nuc = init_state.Tgt().HitNucPdg();
135  int nu = init_state.ProbePdg();
136 
137  bool isP = pdg::IsProton(nuc);
138  bool isN = pdg::IsNeutron(nuc);
139  bool isnu = pdg::IsNeutrino(nu);
140  bool isnub = pdg::IsAntiNeutrino(nu);
141 
142  bool prcok = proc_info.IsWeakCC() && ((isP&&isnub) || (isN&&isnu));
143  if(!prcok) return false;
144 
145  return true;
146 }
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
int HitNucPdg(void) const
Definition: Target.cxx:304
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:338
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:333
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:46
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
int ProbePdg(void) const
Definition: InitialState.h:64
const Target & Tgt(void) const
Definition: InitialState.h:66
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
Initial State information.
Definition: InitialState.h:48
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;
76  // dimension of kine phase space
78  int kpsdim = 1 + std::count(s.begin(), s.end(), ',');
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 
102 
103  // The algorithm computes d^1xsec/dQ2 or d^2xsec/dQ2dv
104  // Check whether variable tranformation is needed
105  if ( kps != kPSQ2fE && kps != kPSQ2vfE )
106  {
107  double J = 1.;
108  if (kpsdim == 1)
110  else if (kpsdim == 2)
112  xsec *= J;
113  }
114 
115  return xsec;
116 
117 }
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
std::string string
Definition: nybbler.cc:12
double d2sQES_dQ2dv_SM(const Interaction *i) const
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
static string AsString(KinePhaseSpace_t kps)
#define pWARN
Definition: Messenger.h:60
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
double dsQES_dQ2_SM(const Interaction *interaction) const
static QCString * s
Definition: config.cpp:1042

Member Data Documentation

double genie::SmithMonizQELCCPXSec::fcosT_k
mutableprivate

Definition at line 96 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fE_lep
mutableprivate

Definition at line 86 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fE_nu
mutableprivate

Definition at line 85 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fF_A
mutableprivate

Definition at line 99 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fF_M
mutableprivate

Definition at line 98 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fF_P
mutableprivate

Definition at line 100 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fF_V
mutableprivate

Definition at line 97 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fFF_A
mutableprivate

Definition at line 103 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fFF_M
mutableprivate

Definition at line 102 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fFF_V
mutableprivate

Definition at line 101 of file SmithMonizQELCCPXSec.h.

QELFormFactors genie::SmithMonizQELCCPXSec::fFormFactors
mutableprivate

Definition at line 78 of file SmithMonizQELCCPXSec.h.

const QELFormFactorsModelI* genie::SmithMonizQELCCPXSec::fFormFactorsModel
private

Definition at line 79 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fk1
mutableprivate

Definition at line 91 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fk2
mutableprivate

Definition at line 92 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fk7
mutableprivate

Definition at line 93 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fm_tar
mutableprivate

Definition at line 89 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fmm_fin
mutableprivate

Definition at line 88 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fmm_ini
mutableprivate

Definition at line 87 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fmm_tar
mutableprivate

Definition at line 90 of file SmithMonizQELCCPXSec.h.

int genie::SmithMonizQELCCPXSec::fn_NT
mutableprivate

Definition at line 82 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fQ2
mutableprivate

Definition at line 83 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fqqv
mutableprivate

Definition at line 95 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fqv
mutableprivate

Definition at line 94 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fv
mutableprivate

Definition at line 84 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fVud2
private

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

Definition at line 81 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fW_1
mutableprivate

Definition at line 104 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fW_2
mutableprivate

Definition at line 105 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fW_3
mutableprivate

Definition at line 106 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fW_4
mutableprivate

Definition at line 107 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fW_5
mutableprivate

Definition at line 108 of file SmithMonizQELCCPXSec.h.

const XSecIntegratorI* genie::SmithMonizQELCCPXSec::fXSecIntegrator
private

Definition at line 80 of file SmithMonizQELCCPXSec.h.

double genie::SmithMonizQELCCPXSec::fXSecScale
private

external xsec scaling factor

Definition at line 77 of file SmithMonizQELCCPXSec.h.

SmithMonizUtils* genie::SmithMonizQELCCPXSec::sm_utils
mutableprivate

Definition at line 70 of file SmithMonizQELCCPXSec.h.


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