Public Member Functions | Private Member Functions | Private Attributes | List of all members
genie::SuSAv2MECPXSec Class Reference

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

#include <SuSAv2MECPXSec.h>

Inheritance diagram for genie::SuSAv2MECPXSec:
genie::XSecAlgorithmI genie::Algorithm

Public Member Functions

 SuSAv2MECPXSec ()
 
 SuSAv2MECPXSec (string config)
 
virtual ~SuSAv2MECPXSec ()
 
double XSec (const Interaction *i, KinePhaseSpace_t k) 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 config)
 
double PairRatio (const Interaction *i) const
 
- 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)
 Load algorithm configuration. More...
 
double Qvalue (const Interaction &interaction) const
 

Private Attributes

double fXSecScale
 External scaling factor for this cross section. More...
 
const genie::HadronTensorModelIfHadronTensorModel
 
string fKFTable
 
double fEbHe
 
double fEbLi
 
double fEbC
 
double fEbO
 
double fEbMg
 
double fEbAr
 
double fEbCa
 
double fEbFe
 
double fEbNi
 
double fEbSn
 
double fEbAu
 
double fEbPb
 
const XSecIntegratorIfXSecIntegrator
 GSL numerical integrator. More...
 
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)
 
- 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 the SuSAv2-MEC model differential cross section. Uses precomputed hadron tensor tables. Is a concrete implementation of the XSecAlgorithmI interface.

Author
Steven Gardiner <gardiner fnal.gov> Fermi National Acclerator Laboratory

Stephen Dolan <stephen.joseph.dolan cern.ch> European Organization for Nuclear Research (CERN)

G::D. Megias et al., "Meson-exchange currents and quasielastic predictions for charged-current neutrino-12C scattering in the superscaling approach," PRD 91 (2015) 073004

November 2, 2018

Copyright (c) 2003-2018, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE

Definition at line 40 of file SuSAv2MECPXSec.h.

Constructor & Destructor Documentation

SuSAv2MECPXSec::SuSAv2MECPXSec ( )

Definition at line 28 of file SuSAv2MECPXSec.cxx.

28  : XSecAlgorithmI("genie::SuSAv2MECPXSec")
29 {
30 }
SuSAv2MECPXSec::SuSAv2MECPXSec ( string  config)

Definition at line 32 of file SuSAv2MECPXSec.cxx.

33  : XSecAlgorithmI("genie::SuSAv2MECPXSec", config)
34 {
35 }
static Config * config
Definition: config.cpp:1054
SuSAv2MECPXSec::~SuSAv2MECPXSec ( )
virtual

Definition at line 37 of file SuSAv2MECPXSec.cxx.

38 {
39 }

Member Function Documentation

void SuSAv2MECPXSec::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 365 of file SuSAv2MECPXSec.cxx.

366 {
367  Algorithm::Configure(config);
368  this->LoadConfig();
369 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void LoadConfig(void)
Load algorithm configuration.
void SuSAv2MECPXSec::Configure ( std::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 371 of file SuSAv2MECPXSec.cxx.

372 {
374  this->LoadConfig();
375 }
static Config * config
Definition: config.cpp:1054
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void LoadConfig(void)
Load algorithm configuration.
double SuSAv2MECPXSec::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 252 of file SuSAv2MECPXSec.cxx.

253 {
254  double xsec = fXSecIntegrator->Integrate(this,interaction);
255  return xsec;
256 }
const XSecIntegratorI * fXSecIntegrator
GSL numerical integrator.
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
void SuSAv2MECPXSec::LoadConfig ( void  )
private

Load algorithm configuration.

Definition at line 377 of file SuSAv2MECPXSec.cxx.

378 {
379  bool good_config = true ;
380  // Cross section scaling factor
381  GetParamDef("MEC-XSecScale", fXSecScale, 1.) ;
382 
383  fHadronTensorModel = dynamic_cast<const HadronTensorModelI*> ( this->SubAlg("HadronTensorAlg") );
384  if( !fHadronTensorModel ) {
385  good_config = false ;
386  LOG("SuSAv2MECPXSec", pERROR) << "The required HadronTensorAlg does not exist. AlgoID is : " << SubAlg("HadronTensorAlg")->Id();
387  }
388 
389  fXSecIntegrator = dynamic_cast<const XSecIntegratorI*> (this->SubAlg("NumericalIntegrationAlg"));
390  if( !fXSecIntegrator ) {
391  good_config = false ;
392  LOG("SuSAv2MECPXSec", pERROR) << "The required NumericalIntegrationAlg does not exist. AlgId is : " << SubAlg("NumericalIntegrationAlg")->Id() ;
393  }
394 
395  //Fermi momentum tables for scaling
396  this->GetParam( "FermiMomentumTable", fKFTable);
397 
398  //binding energy lookups for scaling
399  this->GetParam( "RFG-NucRemovalE@Pdg=1000020040", fEbHe );
400  this->GetParam( "RFG-NucRemovalE@Pdg=1000030060", fEbLi );
401  this->GetParam( "RFG-NucRemovalE@Pdg=1000060120", fEbC );
402  this->GetParam( "RFG-NucRemovalE@Pdg=1000080160", fEbO );
403  this->GetParam( "RFG-NucRemovalE@Pdg=1000120240", fEbMg );
404  this->GetParam( "RFG-NucRemovalE@Pdg=1000180400", fEbAr );
405  this->GetParam( "RFG-NucRemovalE@Pdg=1000200400", fEbCa );
406  this->GetParam( "RFG-NucRemovalE@Pdg=1000260560", fEbFe );
407  this->GetParam( "RFG-NucRemovalE@Pdg=1000280580", fEbNi );
408  this->GetParam( "RFG-NucRemovalE@Pdg=1000501190", fEbSn );
409  this->GetParam( "RFG-NucRemovalE@Pdg=1000791970", fEbAu );
410  this->GetParam( "RFG-NucRemovalE@Pdg=1000822080", fEbPb );
411 
412  // Read optional MECScaleVsW:
413  fMECScaleAlg = nullptr;
414  if( GetConfig().Exists("MECScaleAlg") ) {
415  fMECScaleAlg = dynamic_cast<const XSecScaleI *> ( this->SubAlg("MECScaleAlg") );
416  if( !fMECScaleAlg ) {
417  good_config = false ;
418  LOG("Susav2MECPXSec", pERROR) << "The required MECScaleAlg cannot be casted. AlgID is : " << SubAlg("MECScaleAlg")->Id() ;
419  }
420  }
421 
422  // Read optional QvalueShifter:
423  fQvalueShifter = nullptr;
424  if( GetConfig().Exists("QvalueShifterAlg") ) {
425  fQvalueShifter = dynamic_cast<const QvalueShifter *> ( this->SubAlg("QvalueShifterAlg") );
426  if( !fQvalueShifter ) {
427  good_config = false ;
428  LOG("SuSAv2MECPXSec", pERROR) << "The required QvalueShifterAlg does not exist. AlgId is : " << SubAlg("QvalueShifterAlg")->Id() ;
429  }
430  }
431 
432  if( ! good_config ) {
433  LOG("SuSAv2MECPXSec", pERROR) << "Configuration has failed.";
434  exit(78) ;
435  }
436 
437 }
#define pERROR
Definition: Messenger.h:59
Cross Section Integrator Interface.
const genie::HadronTensorModelI * fHadronTensorModel
This class is responsible to compute a scaling factor for the XSec.
Definition: XSecScaleI.h:25
const XSecIntegratorI * fXSecIntegrator
GSL numerical integrator.
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:97
const QvalueShifter * fQvalueShifter
double fXSecScale
External scaling factor for this cross section.
Creates hadron tensor objects for use in cross section calculations.
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 XSecScaleI * fMECScaleAlg
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345
double SuSAv2MECPXSec::PairRatio ( const Interaction i) const
Todo:
add the different pair configurations for e-scattering

Definition at line 171 of file SuSAv2MECPXSec.cxx.

172 {
173 
174  // Currently we only have the relative pair contributions for C12.
175  // We hope to add mode later, but for the moment assume the relative
176  // contributions are A-independant.
177 
178  int probe_pdg = interaction->InitState().ProbePdg();
179 
180  HadronTensorType_t tensor_type = kHT_Undefined;
181  HadronTensorType_t pn_tensor_type = kHT_Undefined;
182  if ( pdg::IsNeutrino(probe_pdg) || pdg::IsAntiNeutrino(probe_pdg) ) {
183  tensor_type = kHT_MEC_FullAll;
184  pn_tensor_type = kHT_MEC_Fullpn;
185  }
186  else {
187  // If the probe is not a neutrino, assume that it's an electron
188  // For the moment all electron interactions are pp final state
189  tensor_type = kHT_MEC_EM;
190  pn_tensor_type = kHT_MEC_EM;
191  }
192 
193  // The SuSAv2-MEC hadron tensors are defined using the same conventions
194  // as the Valencia MEC model, so we can use the same sort of tensor
195  // object to describe them.
196  const LabFrameHadronTensorI* tensor
198  tensor_type) );
199 
200  const LabFrameHadronTensorI* tensor_pn
202  pn_tensor_type) );
203 
204 
205  /// \todo add the different pair configurations for e-scattering
206 
207  // If retrieving the tensor failed, complain and return zero
208  if ( !tensor ) {
209  LOG("SuSAv2MEC", pWARN) << "Failed to load a hadronic tensor for the"
210  " nuclide " << kPdgTgtC12;
211  return 0.;
212  }
213 
214  // Check that the input kinematical point is within the range
215  // in which hadron tensors are known (for chosen target)
216  double Ev = interaction->InitState().ProbeE(kRfLab);
217  double Tl = interaction->Kine().GetKV(kKVTl);
218  double costl = interaction->Kine().GetKV(kKVctl);
219  double ml = interaction->FSPrimLepton()->Mass();
220  double Q0 = 0.;
221  double Q3 = 0.;
222 
223  genie::utils::mec::Getq0q3FromTlCostl(Tl, costl, Ev, ml, Q0, Q3);
224 
225  double Q0min = tensor->q0Min();
226  double Q0max = tensor->q0Max();
227  double Q3min = tensor->qMagMin();
228  double Q3max = tensor->qMagMax();
229  if (Q0 < Q0min || Q0 > Q0max || Q3 < Q3min || Q3 > Q3max) {
230  return 1.0;
231  }
232 
233  // The Q-Value essentially corrects q0 to account for nuclear
234  // binding energy in the Valencia model but this effect is already
235  // in Guille's tensors so its set it to 0.
236  // However, additional corrections may be necessary:
237  double Delta_Q_value = Qvalue( * interaction ) ;
238 
239  // Compute the cross section using the hadron tensor
240  double xsec_all = tensor->dSigma_dT_dCosTheta_rosenbluth(interaction, Delta_Q_value);
241  double xsec_pn = tensor_pn->dSigma_dT_dCosTheta_rosenbluth(interaction, Delta_Q_value);
242 
243  //hadron tensor precision can sometimes lead to 0 xsec_pn but finite xsec
244  //seems to cause issues downstream ...
245  if(xsec_pn==0) xsec_pn = 0.00001*xsec_all;
246 
247  double ratio = (1e10*xsec_pn)/(1e10*xsec_all);
248 
249  return ratio;
250 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
const genie::HadronTensorModelI * fHadronTensorModel
Abstract interface for an object that computes the elements ( , , etc.) and structure functions ( ...
virtual double qMagMin() const =0
enum genie::HadronTensorType HadronTensorType_t
virtual const HadronTensorI * GetTensor(int tensor_pdg, HadronTensorType_t type) const =0
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual double q0Min() const =0
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
bool Getq0q3FromTlCostl(double Tl, double costl, double Enu, double ml, double &q0, double &q3)
Definition: MECUtils.cxx:121
#define pWARN
Definition: Messenger.h:60
virtual double dSigma_dT_dCosTheta_rosenbluth(const Interaction *interaction, double Q_value) const =0
double Qvalue(const Interaction &interaction) const
const int kPdgTgtC12
Definition: PDGCodes.h:202
virtual double qMagMax() const =0
virtual double q0Max() const =0
double SuSAv2MECPXSec::Qvalue ( const Interaction interaction) const
private
Todo:
Add more hadron tensors so this scaling is not so terrible

Definition at line 281 of file SuSAv2MECPXSec.cxx.

282 {
283  // Get the hadron tensor for the selected nuclide. Check the probe PDG code
284  // to know whether to use the tensor for CC neutrino scattering or for
285  // electron scattering
286  int target_pdg = interaction.InitState().Tgt().Pdg();
287  int probe_pdg = interaction.InitState().ProbePdg();
288  int tensor_pdg = kPdgTgtC12;
289  int A_request = pdg::IonPdgCodeToA(target_pdg);
290 
291  double Eb_tgt=0;
292  double Eb_ten=0;
293 
294  /// \todo Add more hadron tensors so this scaling is not so terrible
295  // At the moment all we have is Carbon so this is all just a place holder ...
296  if ( A_request == 4 ) {
297  Eb_tgt=fEbHe; Eb_ten=fEbC;
298  // This is for helium 4, but use carbon tensor, may not be ideal ...
299  }
300  else if (A_request < 9) {
301  Eb_tgt=fEbLi; Eb_ten=fEbC;
302  }
303  else if (A_request >= 9 && A_request < 15) {
304  Eb_tgt=fEbC; Eb_ten=fEbC;
305  }
306  else if(A_request >= 15 && A_request < 22) {
307  //tensor_pdg = kPdgTgtO16;
308  // Oxygen tensor has some issues - xsec @ 50 GeV = 45.2835 x 1E-38 cm^2
309  // This is ~ 24 times higher than C
310  // I think it's just a missing scale factor but I need to check.
311  Eb_tgt=fEbO; Eb_ten=fEbC;
312  }
313  else if(A_request >= 22 && A_request < 40) {
314  Eb_tgt=fEbMg; Eb_ten=fEbC;
315  }
316  else if(A_request >= 40 && A_request < 56) {
317  Eb_tgt=fEbAr; Eb_ten=fEbC;
318  }
319  else if(A_request >= 56 && A_request < 119) {
320  Eb_tgt=fEbFe; Eb_ten=fEbC;
321  }
322  else if(A_request >= 119 && A_request < 206) {
323  Eb_tgt=fEbSn; Eb_ten=fEbC;
324  }
325  else if(A_request >= 206) {
326  Eb_tgt=fEbPb; Eb_ten=fEbC;
327  }
328 
329  // SD: The Q-Value essentially corrects q0 to account for nuclear
330  // binding energy in the Valencia model but this effect is already
331  // in Guille's tensors so I'll set it to 0.
332  // However, if I want to scale I need to account for the altered
333  // binding energy. To first order I can use the Delta_Q_value for this.
334  // But this is 2p2h - so binding energy counts twice - use 2*1p1h
335  // value (although what should be done here is still not clear).
336 
337  double Delta_Q_value = 2*(Eb_tgt-Eb_ten);
338 
339  // Apply Qvalue relative shift if needed:
340  if( fQvalueShifter ) {
341  // We have the option to add an additional shift on top of the binding energy correction
342  // The QvalueShifter, is a relative shift to the Q_value.
343  // The Q_value was already taken into account in the hadron tensor. Here we recalculate it
344  // to get the right absolute shift.
345  double tensor_Q_value = genie::utils::mec::Qvalue(tensor_pdg,probe_pdg);
346  double total_Q_value = tensor_Q_value + Delta_Q_value ;
347  double Q_value_shift = total_Q_value * fQvalueShifter -> Shift( interaction.InitState().Tgt() ) ;
348  Delta_Q_value += Q_value_shift ;
349  }
350 
351  // We apply an extra Q-value shift here to account for differences between
352  // the 12C EM MEC tensors currently in use (which have a "baked in" Q-value
353  // already incorporated) and the treatment in Guille's thesis. Differences
354  // between the two lead to a few-tens-of-MeV shift in the energy transfer
355  // distribution for EM MEC. The shift is done in terms of the binding energy
356  // value associated with the original tensor (Eb_ten). Corrections for
357  // scaling to a different target are already handled above.
358  // - S. Gardiner, 1 July 2020
359  bool isEM = interaction.ProcInfo().IsEM();
360  if ( isEM ) Delta_Q_value -= 2. * Eb_ten;
361 
362  return Delta_Q_value ;
363 }
int Pdg(void) const
Definition: Target.h:71
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:60
double Qvalue(int targetpdg, int nupdg)
Definition: MECUtils.cxx:164
int ProbePdg(void) const
Definition: InitialState.h:64
bool IsEM(void) const
const QvalueShifter * fQvalueShifter
const int kPdgTgtC12
Definition: PDGCodes.h:202
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
const Target & Tgt(void) const
Definition: InitialState.h:66
bool SuSAv2MECPXSec::ValidProcess ( const Interaction i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 258 of file SuSAv2MECPXSec.cxx.

259 {
260  if ( interaction->TestBit(kISkipProcessChk) ) return true;
261 
262  const ProcessInfo& proc_info = interaction->ProcInfo();
263  if ( !proc_info.IsMEC() ) {
264  return false;
265  }
266 
267  int probe = interaction->InitState().ProbePdg();
268 
269  bool is_nu = pdg::IsNeutrino( probe );
270  bool is_nub = pdg::IsAntiNeutrino( probe );
271  bool is_chgl = pdg::IsChargedLepton( probe );
272 
273  bool prc_ok = ( proc_info.IsWeakCC() && (is_nu || is_nub) )
274  || ( proc_info.IsEM() && is_chgl );
275 
276  if ( !prc_ok ) return false;
277 
278  return true;
279 }
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:98
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
bool IsMEC(void) const
bool IsEM(void) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double SuSAv2MECPXSec::XSec ( const Interaction i,
KinePhaseSpace_t  k 
) const
virtual

Compute the cross section for the input interaction.

Todo:
add the different pair configurations for e-scattering

Implements genie::XSecAlgorithmI.

Definition at line 41 of file SuSAv2MECPXSec.cxx.

43 {
44  // Don't try to do the calculation if we've been handed an interaction that
45  // doesn't make sense
46  if ( !this->ValidProcess(interaction) ) return 0.;
47 
48  // Get the hadron tensor for the selected nuclide. Check the probe PDG code
49  // to know whether to use the tensor for CC neutrino scattering or for
50  // electron scattering
51  int target_pdg = interaction->InitState().Tgt().Pdg();
52  int probe_pdg = interaction->InitState().ProbePdg();
53  int A_request = pdg::IonPdgCodeToA(target_pdg);
54  int Z_request = pdg::IonPdgCodeToZ(target_pdg);
55  bool need_to_scale = false;
56 
57  HadronTensorType_t tensor_type = kHT_Undefined;
58  if ( pdg::IsNeutrino(probe_pdg) || pdg::IsAntiNeutrino(probe_pdg) ) {
59  tensor_type = kHT_MEC_FullAll;
60  //pn_tensor_type = kHT_MEC_Fullpn;
61  //tensor_type = kHT_MEC_FullAll_wImag;
62  //pn_tensor_type = kHT_MEC_FullAll_wImag;
63  }
64  else {
65  // If the probe is not a neutrino, assume that it's an electron
66  // For the moment all electron interactions are pp final state
67  tensor_type = kHT_MEC_EM;
68  //pn_tensor_type = kHT_MEC_EM;
69  }
70 
71  // Currently we only have the relative pair contributions for C12.
72  int tensor_pdg = kPdgTgtC12;
73  if(tensor_pdg != target_pdg) need_to_scale = true;
74 
75  // The SuSAv2-MEC hadron tensors are defined using the same conventions
76  // as the Valencia MEC model, so we can use the same sort of tensor
77  // object to describe them.
78  const LabFrameHadronTensorI* tensor
79  = dynamic_cast<const LabFrameHadronTensorI*>( fHadronTensorModel->GetTensor(tensor_pdg,
80  tensor_type) );
81 
82  /// \todo add the different pair configurations for e-scattering
83 
84  // If retrieving the tensor failed, complain and return zero
85  if ( !tensor ) {
86  LOG("SuSAv2MEC", pWARN) << "Failed to load a hadronic tensor for the"
87  " nuclide " << tensor_pdg;
88  return 0.;
89  }
90 
91  // Check that the input kinematical point is within the range
92  // in which hadron tensors are known (for chosen target)
93  double Ev = interaction->InitState().ProbeE(kRfLab);
94  double Tl = interaction->Kine().GetKV(kKVTl);
95  double costl = interaction->Kine().GetKV(kKVctl);
96  double ml = interaction->FSPrimLepton()->Mass();
97  double Q0 = 0.;
98  double Q3 = 0.;
99 
100  // The Q-Value essentially corrects q0 to account for nuclear
101  // binding energy in the Valencia model but this effect is already
102  // in Guille's tensors so its set it to 0.
103  // However, additional corrections may be necessary:
104  double Delta_Q_value = Qvalue( * interaction ) ;
105 
106  genie::utils::mec::Getq0q3FromTlCostl(Tl, costl, Ev, ml, Q0, Q3);
107 
108  double Q0min = tensor->q0Min();
109  double Q0max = tensor->q0Max();
110  double Q3min = tensor->qMagMin();
111  double Q3max = tensor->qMagMax();
112  if (Q0-Delta_Q_value < Q0min || Q0-Delta_Q_value > Q0max || Q3 < Q3min || Q3 > Q3max) {
113  return 0.0;
114  }
115 
116  // *** Enforce the global Q^2 cut (important for EM scattering) ***
117  // Choose the appropriate minimum Q^2 value based on the interaction
118  // mode (this is important for EM interactions since the differential
119  // cross section blows up as Q^2 --> 0)
120  double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit
121  if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics
122  ::electromagnetic::kMinQ2Limit; // EM limit
123 
124  // Neglect shift due to binding energy. The cut is on the actual
125  // value of Q^2, not the effective one to use in the tensor contraction.
126  double Q2 = Q3*Q3 - Q0*Q0;
127  if ( Q2 < Q2min ) return 0.;
128 
129  // By default, we will compute the full cross-section. If a {p,n} hit
130  // dinucleon was set we will calculate the cross-section for that
131  // component only
132 
133  bool pn = (interaction->InitState().Tgt().HitNucPdg() == kPdgClusterNP);
134 
135  // Compute the cross section using the hadron tensor
136  double xsec = tensor->dSigma_dT_dCosTheta_rosenbluth(interaction, Delta_Q_value);
137 
138  // This scaling should be okay-ish for the total xsec, but it misses
139  // the energy shift. To get this we should really just build releveant
140  // hadron tensors but there may be some ways to approximate it.
141  // For more details see Guille's thesis: https://idus.us.es/xmlui/handle/11441/74826
142  if ( need_to_scale ) {
144  const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
145  double KF_tgt = kft->FindClosestKF(target_pdg, kPdgProton);
146  double KF_ten = kft->FindClosestKF(tensor_pdg, kPdgProton);
147  LOG("SuSAv2MEC", pDEBUG) << "KF_tgt = " << KF_tgt;
148  LOG("SuSAv2MEC", pDEBUG) << "KF_ten = " << KF_ten;
149  double A_ten = pdg::IonPdgCodeToA(tensor_pdg);
150  double scaleFact = (A_request/A_ten)*(KF_tgt/KF_ten)*(KF_tgt/KF_ten);
151  xsec *= scaleFact;
152  }
153 
154  // Apply given overall scaling factor
155  xsec *= fXSecScale;
156 
157  // Scale given a scaling algorithm:
158  if( fMECScaleAlg ) xsec *= fMECScaleAlg->GetScaling( * interaction ) ;
159 
160  if ( kps != kPSTlctl ) {
161  LOG("SuSAv2MEC", pWARN)
162  << "Doesn't support transformation from "
163  << KinePhaseSpace::AsString(kPSTlctl) << " to "
164  << KinePhaseSpace::AsString(kps);
165  xsec = 0.;
166  }
167 
168  return xsec;
169 }
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:107
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1064
const genie::HadronTensorModelI * fHadronTensorModel
const int kPdgClusterNP
Definition: PDGCodes.h:215
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:60
static FermiMomentumTablePool * Instance(void)
Abstract interface for an object that computes the elements ( , , etc.) and structure functions ( ...
A table of Fermi momentum constants.
virtual double qMagMin() const =0
enum genie::HadronTensorType HadronTensorType_t
virtual const HadronTensorI * GetTensor(int tensor_pdg, HadronTensorType_t type) const =0
static const double kMinQ2Limit
Definition: Controls.h:41
Singleton class to load & serve tables of Fermi momentum constants.
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
virtual double q0Min() const =0
const FermiMomentumTable * GetTable(string name)
virtual double GetScaling(const Interaction &) const =0
static string AsString(KinePhaseSpace_t kps)
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:115
bool Getq0q3FromTlCostl(double Tl, double costl, double Enu, double ml, double &q0, double &q3)
Definition: MECUtils.cxx:121
#define pWARN
Definition: Messenger.h:60
virtual double dSigma_dT_dCosTheta_rosenbluth(const Interaction *interaction, double Q_value) const =0
double Qvalue(const Interaction &interaction) const
Kinematical utilities.
const int kPdgTgtC12
Definition: PDGCodes.h:202
double fXSecScale
External scaling factor for this cross section.
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
virtual double qMagMax() const =0
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:52
const int kPdgProton
Definition: PDGCodes.h:81
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
const XSecScaleI * fMECScaleAlg
virtual double q0Max() const =0
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

double genie::SuSAv2MECPXSec::fEbAr
private

Definition at line 84 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbAu
private

Definition at line 89 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbC
private

Definition at line 81 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbCa
private

Definition at line 85 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbFe
private

Definition at line 86 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbHe
private

Definition at line 79 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbLi
private

Definition at line 80 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbMg
private

Definition at line 83 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbNi
private

Definition at line 87 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbO
private

Definition at line 82 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbPb
private

Definition at line 90 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fEbSn
private

Definition at line 88 of file SuSAv2MECPXSec.h.

const genie::HadronTensorModelI* genie::SuSAv2MECPXSec::fHadronTensorModel
private

Definition at line 73 of file SuSAv2MECPXSec.h.

string genie::SuSAv2MECPXSec::fKFTable
private

Definition at line 76 of file SuSAv2MECPXSec.h.

const XSecScaleI* genie::SuSAv2MECPXSec::fMECScaleAlg
private

Definition at line 95 of file SuSAv2MECPXSec.h.

const QvalueShifter* genie::SuSAv2MECPXSec::fQvalueShifter
private

Definition at line 96 of file SuSAv2MECPXSec.h.

const XSecIntegratorI* genie::SuSAv2MECPXSec::fXSecIntegrator
private

GSL numerical integrator.

Definition at line 93 of file SuSAv2MECPXSec.h.

double genie::SuSAv2MECPXSec::fXSecScale
private

External scaling factor for this cross section.

Definition at line 71 of file SuSAv2MECPXSec.h.


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