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

Computes the Valencia MEC model differential cross section. Uses precomputed hadon tensor tables. Is a concrete implementation of the XSecAlgorithmI interface. More...

#include <NievesSimoVacasMECPXSec2016.h>

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

Public Member Functions

 NievesSimoVacasMECPXSec2016 ()
 NievesSimoVacasMECPXSec2016 (string config)
virtual ~NievesSimoVacasMECPXSec2016 ()
double XSec (const Interaction *i, KinePhaseSpace_t k) 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 config)
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)

Private Attributes

double fXSecCCScale
 external xsec scaling factor
double fXSecNCScale
 external xsec scaling factor
const HadronTensorModelIfHadronTensorModel
const XSecIntegratorIfXSecIntegrator
const XSecScaleIfMECScaleAlg
const QvalueShifterfQvalueShifter

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 the Valencia MEC model differential cross section. Uses precomputed hadon tensor tables. Is a concrete implementation of the XSecAlgorithmI interface.

Author
Code contributed by J. Schwehr, D. Cherdack, R. Gran and described in arXiv:1601.02038 and some of the refereces there-in, in particular PRD 88 (2013) 113007

Substantial code refactorizations by the core GENIE group.

Refactored in 2018 by S. Gardiner to use the new hadron tensor framework

References:\n J. Nieves, I. Ruiz Simo, M.J. Vicente Vacas,
Inclusive quasi-elastic neutrino reactions, PRC 83 (2011) 045501
Created:\n Mar 22, 2016
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 41 of file NievesSimoVacasMECPXSec2016.h.

Constructor & Destructor Documentation

◆ NievesSimoVacasMECPXSec2016() [1/2]

NievesSimoVacasMECPXSec2016::NievesSimoVacasMECPXSec2016 ( )

Definition at line 29 of file NievesSimoVacasMECPXSec2016.cxx.

29 :
30XSecAlgorithmI("genie::NievesSimoVacasMECPXSec2016")
31{
32
33}

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ NievesSimoVacasMECPXSec2016() [2/2]

NievesSimoVacasMECPXSec2016::NievesSimoVacasMECPXSec2016 ( string config)

Definition at line 35 of file NievesSimoVacasMECPXSec2016.cxx.

35 :
36XSecAlgorithmI("genie::NievesSimoVacasMECPXSec2016", config)
37{
38
39}

References genie::XSecAlgorithmI::XSecAlgorithmI().

◆ ~NievesSimoVacasMECPXSec2016()

NievesSimoVacasMECPXSec2016::~NievesSimoVacasMECPXSec2016 ( )
virtual

Definition at line 41 of file NievesSimoVacasMECPXSec2016.cxx.

42{
43
44}

Member Function Documentation

◆ Configure() [1/2]

void NievesSimoVacasMECPXSec2016::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 360 of file NievesSimoVacasMECPXSec2016.cxx.

361{
362 Algorithm::Configure(config);
363 this->LoadConfig();
364}
virtual void Configure(const Registry &config)
Definition Algorithm.cxx:62

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

◆ Configure() [2/2]

void NievesSimoVacasMECPXSec2016::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 366 of file NievesSimoVacasMECPXSec2016.cxx.

367{
368 Algorithm::Configure(config);
369 this->LoadConfig();
370}

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

◆ Integral()

double NievesSimoVacasMECPXSec2016::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 341 of file NievesSimoVacasMECPXSec2016.cxx.

343{
344 double xsec = fXSecIntegrator->Integrate(this,interaction);
345 return xsec;
346}

References fXSecIntegrator.

◆ LoadConfig()

void NievesSimoVacasMECPXSec2016::LoadConfig ( void )
private

Definition at line 372 of file NievesSimoVacasMECPXSec2016.cxx.

373{
374 bool good_config = true;
375
376 // Cross section scaling factor
377 GetParam( "MEC-CC-XSecScale", fXSecCCScale ) ;
378 GetParam( "MEC-NC-XSecScale", fXSecNCScale ) ;
379
380 fHadronTensorModel = dynamic_cast<const HadronTensorModelI *> ( this->SubAlg("HadronTensorAlg") );
381 if( !fHadronTensorModel ) {
382 good_config = false ;
383 LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required HadronTensorAlg does not exist. AlgID is : " << SubAlg("HadronTensorAlg")->Id() ;
384 }
385
386 fXSecIntegrator = dynamic_cast<const XSecIntegratorI *> (this->SubAlg("NumericalIntegrationAlg"));
387 if( !fXSecIntegrator ) {
388 good_config = false ;
389 LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required NumericalIntegrationAlg does not exist. AlgID is : " << SubAlg("NumericalIntegrationAlg")->Id();
390 }
391
392 // Read optional QvalueShifter:
393 fQvalueShifter = nullptr;
394 if( GetConfig().Exists("QvalueShifterAlg") ) {
395 fQvalueShifter = dynamic_cast<const QvalueShifter *> ( this->SubAlg("QvalueShifterAlg") );
396 if( !fQvalueShifter ) {
397 good_config = false ;
398 LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required QvalueShifterAlg does not exist. AlgID is : " << SubAlg("QvalueShifterAlg")->Id() ;
399 }
400 }
401
402 // Read optional MECScaleVsW:
403 fMECScaleAlg = nullptr;
404 if( GetConfig().Exists("MECScaleAlg") ) {
405 fMECScaleAlg = dynamic_cast<const XSecScaleI *> ( this->SubAlg("MECScaleAlg") );
406 if( !fMECScaleAlg ) {
407 good_config = false ;
408 LOG("NievesSimoVacasMECPXSec2016", pERROR) << "The required MECScaleAlg cannot be casted. AlgID is : " << SubAlg("MECScaleAlg")->Id() ;
409 }
410 }
411
412 if( ! good_config ) {
413 LOG("NievesSimoVacasMECPXSec2016", pERROR) << "Configuration has failed.";
414 exit(78) ;
415 }
416
417}
#define pERROR
Definition Messenger.h:59
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE,...
Definition Messenger.h:96
virtual const Registry & GetConfig(void) const
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Algorithm * SubAlg(const RgKey &registry_key) const
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition Algorithm.h:98
double fXSecCCScale
external xsec scaling factor
double fXSecNCScale
external xsec scaling factor

References fHadronTensorModel, fMECScaleAlg, fQvalueShifter, fXSecCCScale, fXSecIntegrator, fXSecNCScale, genie::Algorithm::GetConfig(), genie::Algorithm::GetParam(), genie::Algorithm::Id(), LOG, pERROR, and genie::Algorithm::SubAlg().

Referenced by Configure(), and Configure().

◆ ValidProcess()

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

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 348 of file NievesSimoVacasMECPXSec2016.cxx.

350{
351 if (interaction->TestBit(kISkipProcessChk)) return true;
352
353 const ProcessInfo & proc_info = interaction->ProcInfo();
354 if (!proc_info.IsMEC()) {
355 return false;
356 }
357 return true;
358}
bool IsMEC(void) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition Interaction.h:47

References genie::ProcessInfo::IsMEC(), genie::kISkipProcessChk, and genie::Interaction::ProcInfo().

◆ XSec()

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

Compute the cross section for the input interaction.

Todo
Replace these hard-coded replacements with an equivalent XML configuration
Todo
Revisit this assumption, and perhaps implement something more robust
Todo
Shouldn't we get this from the nuclear model?

Implements genie::XSecAlgorithmI.

Definition at line 46 of file NievesSimoVacasMECPXSec2016.cxx.

48{
49 // If {W,Q2} have been supplied instead, compute {Tl, ctl}
50 // NOTE: The expressions used here neglect Fermi motion and
51 // should eventually be revisited. See the "important note"
52 // in src/Framework/Utils/KineUtils.cxx about the
53 // Jacobian for transforming {W,Q2} --> {Tl, ctl}.
54 // - S. Gardiner, 29 July 2020
55 if ( kps == kPSWQ2fE ) {
56
57 double Q2 = interaction->Kine().GetKV( kKVQ2 );
58 double W = interaction->Kine().GetKV( kKVW );
59
60 // Probe properties (mass, energy, momentum)
61 const InitialState& init_state = interaction->InitState();
62 double mv = init_state.Probe()->Mass();
63 double Ev = init_state.ProbeE( kRfLab );
64 double pv = std::sqrt( std::max(0., Ev*Ev - mv*mv) );
65
66 // Invariant mass of the initial hit nucleon
67 const TLorentzVector& hit_nuc_P4 = init_state.Tgt().HitNucP4();
68 double M = hit_nuc_P4.M();
69
70 // Get the outgoing lepton kinetic energy
71 double ml = interaction->FSPrimLepton()->Mass();
72 double Tl = Ev - ml - ( (W*W + Q2 - M*M) / (2.*M) );
73
74 // Get the outgoing lepton scattering cosine
75 double El = Tl + ml;
76 double pl = std::sqrt( std::max(0., El*El - ml*ml) );
77 double ctl = ( 2.*Ev*El - Q2 - mv*mv - ml*ml ) / ( 2. * pv * pl );
78
79 // Set Tl, ctl in the interaction
80 interaction->KinePtr()->SetKV( kKVTl, Tl );
81 interaction->KinePtr()->SetKV( kKVctl, ctl );
82 }
83
84 // This function returns d2sigma/(dTmu dcos_mu) in GeV^(-3)
85 int target_pdg = interaction->InitState().Tgt().Pdg();
86
87 int A_request = pdg::IonPdgCodeToA(target_pdg);
88 int Z_request = pdg::IonPdgCodeToZ(target_pdg);
89
90 // To generate cross-sections for nuclei other than those with hadron
91 // tensors we need to pull both the full cross-section and
92 // the pn initial state fraction.
93 // Non-isoscalar nuclei are beyond the original published Valencia model
94 // and scale with A according to the number of pp, pn, or nn pairs
95 // the probe is expected to find.
96 // There is some by-hand optimization here, skipping the delta part when
97 // only the total cross-section is requested.
98 // Possible future models without a Delta had tensor would also use that
99 // flag to call this without computing the Delta part.
100
101 // Try to look up a hadron tensor in the pool that is an exact match for
102 // the target nucleus. If an exact match cannot be found, decide upon a
103 // suitable substitute based on the mass number A and proton number Z.
104
105 int tensor_pdg = target_pdg;
106
107 /// \todo Replace these hard-coded replacements with an equivalent XML
108 /// configuration
109 if ( ! fHadronTensorModel->GetTensor(tensor_pdg, genie::kHT_MEC_FullAll) )
110 {
111
112 if ( A_request == 4 && Z_request == 2 ) {
113 tensor_pdg = kPdgTgtC12;
114 // This is for helium 4, but use carbon tensor
115 // the use of nuclear density parameterization is suspicious
116 // but some users (MINERvA) need something not nothing.
117 // The pn will be exactly 1/3, but pp and nn will be ~1/4
118 // Because the combinatorics are different.
119 // Could do lithium beryllium boron which you don't need
120 }
121 else if (A_request < 9) {
122 // refuse to do D, T, He3, Li, and some Be, B
123 // actually it would work technically, maybe except D, T
124 MAXLOG("NievesSimoVacasMEC", pWARN, 10)
125 << "Asked to scale to deuterium through boron "
126 << target_pdg << " nope, lets not do that.";
127 return 0;
128 }
129 else if (A_request >= 9 && A_request < 15) {
130 tensor_pdg = kPdgTgtC12;
131 //}
132 // could explicitly put in nitrogen for air
133 //else if ( A_request >= 14 && A < 15) { // AND CHANGE <=14 to <14.
134 // tensor_pdg = kPdgTgtN14;
135 }
136 else if (A_request >= 15 && A_request < 22) {
137 tensor_pdg = kPdgTgtO16;
138 }
139 else if (A_request >= 22 && A_request < 33) {
140 // of special interest, this gets Al27 and Si28
141 tensor_pdg = 1000140280;
142 }
143 else if (A_request >= 33 && A_request < 50) {
144 // of special interest, this gets Ar40 and Ti48
145 tensor_pdg = kPdgTgtCa40;
146 }
147 else if (A_request >= 50 && A_request < 90) {
148 // pseudoFe56, also covers many other ferrometals and Ge
149 tensor_pdg = 1000280560;
150 }
151 else if (A_request >= 90 && A_request < 160) {
152 // use Ba112 = PseudoCd. Row5 of Periodic table useless. Ag, Xe?
153 tensor_pdg = 1000561120;
154 }
155 else if (A_request >= 160) {
156 // use Rf208 = pseudoPb
157 tensor_pdg = 1001042080;
158 }
159 else {
160 MAXLOG("NievesSimoVacasMEC", pWARN, 10)
161 << "Asked to scale to a nucleus "
162 << target_pdg << " which we don't know yet.";
163 return 0;
164 }
165 }
166
167 // Check that the input kinematical point is within the range
168 // in which hadron tensors are known (for chosen target)
169 double Ev = interaction->InitState().ProbeE(kRfLab);
170 double Tl = interaction->Kine().GetKV(kKVTl);
171 double costl = interaction->Kine().GetKV(kKVctl);
172 double ml = interaction->FSPrimLepton()->Mass();
173 double Q0 = interaction->Kine().GetKV(kKVQ0);
174 double Q3 = interaction->Kine().GetKV(kKVQ3);
175
176 const LabFrameHadronTensorI* tensor
177 = dynamic_cast<const LabFrameHadronTensorI*>(
178 fHadronTensorModel->GetTensor(tensor_pdg, genie::kHT_MEC_FullAll) );
179
180 // If retrieving the tensor failed, complain and return zero
181 if ( !tensor ) {
182 LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a"
183 " hadronic tensor for the nuclide " << tensor_pdg;
184 return 0.;
185 }
186
187 // Assume for now that the range of validity for the "FullAll" hadron
188 // tensor is the same as for the partial hadron tensors
189 /// \todo Revisit this assumption, and perhaps implement something
190 /// more robust
191 double Q0min = tensor->q0Min();
192 double Q0max = tensor->q0Max();
193 double Q3min = tensor->qMagMin();
194 double Q3max = tensor->qMagMax();
195 if (Q0 < Q0min || Q0 > Q0max || Q3 < Q3min || Q3 > Q3max) {
196 return 0.0;
197 }
198
199 // Get the Q-value needed to calculate the cross sections using the
200 // hadron tensor.
201 /// \todo Shouldn't we get this from the nuclear model?
202 int nu_pdg = interaction->InitState().ProbePdg();
203 double Q_value = genie::utils::mec::Qvalue(target_pdg, nu_pdg);
204
205 // Apply Qvalue relative shift if needed:
206 if( fQvalueShifter ) Q_value += Q_value * fQvalueShifter -> Shift( interaction->InitState().Tgt() ) ;
207
208 // By default, we will compute the full cross-section. If a resonance is
209 // set, we will calculate the part of the cross-section with an internal
210 // Delta line without a final state pion (usually called PPD for pioness
211 // Delta decay). If a {p,n} hit dinucleon was set we will calculate the
212 // cross-section for that component only (either full or PDD cross-section)
213 bool delta = interaction->ExclTag().KnownResonance();
214 bool pn = (interaction->InitState().Tgt().HitNucPdg() == kPdgClusterNP);
215
216 double xsec_all = 0.;
217 double xsec_pn = 0.;
218 if ( delta ) {
219
220 const LabFrameHadronTensorI* tensor_delta_all
221 = dynamic_cast<const LabFrameHadronTensorI*>(
222 fHadronTensorModel->GetTensor(tensor_pdg, genie::kHT_MEC_DeltaAll) );
223
224 if ( !tensor_delta_all ) {
225 LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"DeltaAll\""
226 << " hadronic tensor for nuclide " << tensor_pdg;
227 return 0.;
228 }
229
230 const LabFrameHadronTensorI* tensor_delta_pn
231 = dynamic_cast<const LabFrameHadronTensorI*>(
232 fHadronTensorModel->GetTensor(tensor_pdg, genie::kHT_MEC_Deltapn) );
233
234 if ( !tensor_delta_pn ) {
235 LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"Deltapn\""
236 << " hadronic tensor for nuclide " << tensor_pdg;
237 return 0.;
238 }
239
240 xsec_all = tensor_delta_all->dSigma_dT_dCosTheta(interaction, Q_value);
241 xsec_pn = tensor_delta_pn->dSigma_dT_dCosTheta(interaction, Q_value);
242
243 }
244 else {
245
246 const LabFrameHadronTensorI* tensor_full_all
247 = dynamic_cast<const LabFrameHadronTensorI*>(
248 fHadronTensorModel->GetTensor(tensor_pdg, genie::kHT_MEC_FullAll) );
249
250 if ( !tensor_full_all ) {
251 LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"FullAll\""
252 << " hadronic tensor for nuclide " << tensor_pdg;
253 return 0.;
254 }
255
256 const LabFrameHadronTensorI* tensor_full_pn
257 = dynamic_cast<const LabFrameHadronTensorI*>(
258 fHadronTensorModel->GetTensor(tensor_pdg, genie::kHT_MEC_Fullpn) );
259
260 if ( !tensor_full_pn ) {
261 LOG("NievesSimoVacasMEC", pWARN) << "Failed to load a \"Fullpn\""
262 << " hadronic tensor for nuclide " << tensor_pdg;
263 return 0.;
264 }
265
266 xsec_all = tensor_full_all->dSigma_dT_dCosTheta(interaction, Q_value);
267 xsec_pn = tensor_full_pn->dSigma_dT_dCosTheta(interaction, Q_value);
268 }
269
270 // We need to scale the cross section appropriately if
271 // we are using a hadronic tensor for a nuclide that is different
272 // from the actual target
273 bool need_to_scale = (target_pdg != tensor_pdg);
274
275 // would need to trap and treat He3, T, D special here.
276 if ( need_to_scale ) {
277
278 double PP = Z_request;
279 double NN = A_request - PP;
280 double P = pdg::IonPdgCodeToZ(tensor_pdg);
281 double N = pdg::IonPdgCodeToA(tensor_pdg) - P;
282
283 double scale_pn = TMath::Sqrt( (PP*NN)/(P*N) );
284 double scale_pp = TMath::Sqrt( (PP * (PP - 1.)) / (P * (P - 1.)) );
285 double scale_nn = TMath::Sqrt( (NN * (NN - 1.)) / (N * (N - 1.)) );
286
287 LOG("NievesSimoVacasMEC", pDEBUG)
288 << "Scale pn pp nn for (" << target_pdg << ", " << tensor_pdg << ")"
289 << " : " << scale_pn << " " << scale_pp << " " << scale_nn;
290
291 // This is an approximation in at least three senses:
292 // 1. We are scaling from an isoscalar nucleus using p and n counting
293 // 2. We are not using the right qvalue in the had tensor
294 // 3. We are not scaling the Delta faster than the non-Delta.
295 // The guess is that these are good approximations.
296 // A test we could document is to scale from O16 to N14 or C12 using this
297 // algorithm and see how many percent deviation we see from the full
298 // calculation.
299 double temp_all = xsec_all;
300 double temp_pn = xsec_pn * scale_pn;
301 if (nu_pdg > 0) {
302 // matter neutrinos
303 temp_all = xsec_pn * scale_pn + (xsec_all - xsec_pn) * scale_nn;
304 }
305 else {
306 // antineutrinos
307 temp_all = xsec_pn * scale_pn + (xsec_all - xsec_pn) * scale_pp;
308 }
309
310 xsec_all = temp_all;
311 xsec_pn = temp_pn;
312
313 }
314
315 // Choose the right kind of cross section ("all" or "pn") to return
316 // based on whether a {p, n} dinucleon was hit
317 double xsec = (pn) ? xsec_pn : xsec_all;
318
319 // Apply given scaling factor
320 const ProcessInfo& proc_info = interaction->ProcInfo();
321 if( proc_info.IsWeakCC() ) xsec *= fXSecCCScale;
322 else if( proc_info.IsWeakNC() ) xsec *= fXSecNCScale;
323
324 if( fMECScaleAlg ) xsec *= fMECScaleAlg->GetScaling( * interaction ) ;
325
326 if ( kps != kPSTlctl && kps != kPSWQ2fE ) {
327 LOG("NievesSimoVacasMEC", pWARN)
328 << "Doesn't support transformation from "
331 xsec = 0;
332 }
333 else if ( kps == kPSWQ2fE && xsec != 0. ) {
334 double J = utils::kinematics::Jacobian( interaction, kPSTlctl, kps );
335 xsec *= J;
336 }
337
338 return xsec;
339}
#define MAXLOG(s, l, c)
Similar to LOG(stream,priority) but quits after "maxcount" messages.
Definition Messenger.h:241
#define pDEBUG
Definition Messenger.h:63
#define pWARN
Definition Messenger.h:60
virtual double q0Min() const =0
virtual double qMagMax() const =0
virtual double q0Max() const =0
virtual double qMagMin() const =0
const Target & Tgt(void) const
TParticlePDG * Probe(void) const
double ProbeE(RefFrame_t rf) const
static string AsString(KinePhaseSpace_t kps)
virtual double dSigma_dT_dCosTheta(const Interaction *interaction, double Q_value) const =0
bool IsWeakNC(void) const
bool IsWeakCC(void) const
const TLorentzVector & HitNucP4(void) const
Definition Target.h:91
int IonPdgCodeToZ(int pdgc)
Definition PDGUtils.cxx:55
int IonPdgCodeToA(int pdgc)
Definition PDGUtils.cxx:63
double W(const Interaction *const i)
double Q2(const Interaction *const i)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
double Qvalue(int targetpdg, int nupdg)
Definition MECUtils.cxx:164
double J(double q0, double q3, double Enu, double ml)
Definition MECUtils.cxx:147
const int kPdgClusterNP
Definition PDGCodes.h:215
@ kHT_MEC_Deltapn
@ kHT_MEC_FullAll
@ kHT_MEC_DeltaAll
@ kHT_MEC_Fullpn
const int kPdgTgtCa40
Definition PDGCodes.h:204
@ kKVQ3
Definition KineVar.h:58
@ kKVQ2
Definition KineVar.h:33
@ kKVTl
Definition KineVar.h:38
@ kKVW
Definition KineVar.h:35
@ kKVctl
Definition KineVar.h:39
@ kKVQ0
Definition KineVar.h:57
@ kRfLab
Definition RefFrame.h:26
const int kPdgTgtO16
Definition PDGCodes.h:203
const int kPdgTgtC12
Definition PDGCodes.h:202

References genie::KinePhaseSpace::AsString(), genie::LabFrameHadronTensorI::dSigma_dT_dCosTheta(), genie::Interaction::ExclTag(), fHadronTensorModel, fMECScaleAlg, fQvalueShifter, genie::Interaction::FSPrimLepton(), fXSecCCScale, fXSecNCScale, genie::Kinematics::GetKV(), genie::Target::HitNucP4(), genie::Interaction::InitState(), genie::pdg::IonPdgCodeToA(), genie::pdg::IonPdgCodeToZ(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), genie::utils::kinematics::Jacobian(), genie::kHT_MEC_DeltaAll, genie::kHT_MEC_Deltapn, genie::kHT_MEC_FullAll, genie::kHT_MEC_Fullpn, genie::Interaction::Kine(), genie::Interaction::KinePtr(), genie::kKVctl, genie::kKVQ0, genie::kKVQ2, genie::kKVQ3, genie::kKVTl, genie::kKVW, genie::XclsTag::KnownResonance(), genie::kPdgClusterNP, genie::kPdgTgtC12, genie::kPdgTgtCa40, genie::kPdgTgtO16, genie::kPSTlctl, genie::kPSWQ2fE, genie::kRfLab, LOG, MAXLOG, pDEBUG, genie::Target::Pdg(), genie::InitialState::Probe(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), pWARN, genie::HadronTensorI::q0Max(), genie::HadronTensorI::q0Min(), genie::HadronTensorI::qMagMax(), genie::HadronTensorI::qMagMin(), genie::utils::mec::Qvalue(), genie::Kinematics::SetKV(), and genie::InitialState::Tgt().

Member Data Documentation

◆ fHadronTensorModel

const HadronTensorModelI* genie::NievesSimoVacasMECPXSec2016::fHadronTensorModel
private

Definition at line 66 of file NievesSimoVacasMECPXSec2016.h.

Referenced by LoadConfig(), and XSec().

◆ fMECScaleAlg

const XSecScaleI* genie::NievesSimoVacasMECPXSec2016::fMECScaleAlg
private

Definition at line 70 of file NievesSimoVacasMECPXSec2016.h.

Referenced by LoadConfig(), and XSec().

◆ fQvalueShifter

const QvalueShifter* genie::NievesSimoVacasMECPXSec2016::fQvalueShifter
private

Definition at line 71 of file NievesSimoVacasMECPXSec2016.h.

Referenced by LoadConfig(), and XSec().

◆ fXSecCCScale

double genie::NievesSimoVacasMECPXSec2016::fXSecCCScale
private

external xsec scaling factor

Definition at line 63 of file NievesSimoVacasMECPXSec2016.h.

Referenced by LoadConfig(), and XSec().

◆ fXSecIntegrator

const XSecIntegratorI* genie::NievesSimoVacasMECPXSec2016::fXSecIntegrator
private

Definition at line 68 of file NievesSimoVacasMECPXSec2016.h.

Referenced by Integral(), and LoadConfig().

◆ fXSecNCScale

double genie::NievesSimoVacasMECPXSec2016::fXSecNCScale
private

external xsec scaling factor

Definition at line 64 of file NievesSimoVacasMECPXSec2016.h.

Referenced by LoadConfig(), and XSec().


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