GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
Public Member Functions | Private Types | Private Member Functions | Private Attributes | List of all members
genie::SuSAv2QELPXSec Class Reference

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

#include <SuSAv2QELPXSec.h>

Inheritance diagram for genie::SuSAv2QELPXSec:
Inheritance graph
[legend]
Collaboration diagram for genie::SuSAv2QELPXSec:
Collaboration graph
[legend]

Public Member Functions

 SuSAv2QELPXSec ()
 
 SuSAv2QELPXSec (string config)
 
virtual ~SuSAv2QELPXSec ()
 
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)
 
- 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 Types

enum  modelType {
  kMd_Undefined = 0, kMd_SuSAv2 = 1, kMd_CRPA = 2, kMd_HF = 3,
  kMd_CRPASuSAv2Hybrid = 4, kMd_HFSuSAv2Hybrid = 5, kMd_CRPAPW = 6, kMd_HFPW = 7,
  kMd_CRPAPWSuSAv2Hybrid = 8, kMd_HFPWSuSAv2Hybrid = 9, kMd_SuSAv2Blend = 10
}
 

Private Member Functions

void LoadConfig (void)
 Load algorithm configuration. More...
 
double XSecScaling (double xsec, const Interaction *i, int target_pdg, int tensor_pdg, bool need_to_scale) const
 

Private Attributes

double fXSecCCScale
 External scaling factor for this cross section. More...
 
double fXSecNCScale
 
double fXSecEMScale
 
int blendMode
 Blending start/end q0 value for the combined models. More...
 
double q0BlendStart
 
double q0BlendEnd
 
double qBlendDel
 
double qBlendRef
 
modelType modelConfig
 
const 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 XSecAlgorithmIfFreeNucleonXSecAlg
 Alternate cross section model for free nucleon targets. More...
 
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. 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
 
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. More...
 
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 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< bool > fOwnerships
 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-QE 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)

References:
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
Created:
November 2, 2018
License:
Copyright (c) 2003-2024, 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 SuSAv2QELPXSec.h.

Member Enumeration Documentation

Constructor & Destructor Documentation

SuSAv2QELPXSec::SuSAv2QELPXSec ( )

Definition at line 33 of file SuSAv2QELPXSec.cxx.

33  : XSecAlgorithmI("genie::SuSAv2QELPXSec")
34 {
35 }
SuSAv2QELPXSec::SuSAv2QELPXSec ( string  config)

Definition at line 37 of file SuSAv2QELPXSec.cxx.

38  : XSecAlgorithmI("genie::SuSAv2QELPXSec", config)
39 {
40 }
SuSAv2QELPXSec::~SuSAv2QELPXSec ( )
virtual

Definition at line 42 of file SuSAv2QELPXSec.cxx.

43 {
44 }

Member Function Documentation

void SuSAv2QELPXSec::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 588 of file SuSAv2QELPXSec.cxx.

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

589 {
590  Algorithm::Configure(config);
591  this->LoadConfig();
592 }
void LoadConfig(void)
Load algorithm configuration.
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void genie::SuSAv2QELPXSec::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.

double SuSAv2QELPXSec::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 553 of file SuSAv2QELPXSec.cxx.

References fXSecIntegrator, and genie::XSecIntegratorI::Integrate().

554 {
555  double xsec = fXSecIntegrator->Integrate(this, interaction);
556  return xsec;
557 }
const XSecIntegratorI * fXSecIntegrator
GSL numerical integrator.
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
void SuSAv2QELPXSec::LoadConfig ( void  )
private

Load algorithm configuration.

Definition at line 600 of file SuSAv2QELPXSec.cxx.

References blendMode, fEbAr, fEbAu, fEbC, fEbCa, fEbFe, fEbHe, fEbLi, fEbMg, fEbNi, fEbO, fEbPb, fEbSn, fHadronTensorModel, fKFTable, fQvalueShifter, fXSecCCScale, fXSecEMScale, fXSecIntegrator, fXSecNCScale, genie::Algorithm::GetConfig(), genie::Algorithm::GetParam(), genie::Algorithm::Id(), LOG, modelConfig, pERROR, q0BlendEnd, q0BlendStart, qBlendDel, qBlendRef, and genie::Algorithm::SubAlg().

Referenced by Configure().

601 {
602  bool good_config = true ;
603 
604  // Cross section scaling factor
605  GetParam( "QEL-CC-XSecScale", fXSecCCScale ) ;
606  GetParam( "QEL-NC-XSecScale", fXSecNCScale ) ;
607  GetParam( "QEL-EM-XSecScale", fXSecEMScale ) ;
608 
609  // Cross section model choice
610  int modelChoice;
611  GetParam( "Model-Config", modelChoice ) ;
612  modelConfig = (modelType)modelChoice;
613 
614  // Blending parameters
615  GetParam( "Blend-Mode", blendMode ) ; // 1 = linear, 2 = exp
616  GetParam( "q0-Blend-Start", q0BlendStart ) ; // Used for linear
617  GetParam( "q0-Blend-End", q0BlendEnd ) ; // Used for linear
618  GetParam( "q-Blend-del", qBlendDel ) ; // Used for exp
619  GetParam( "q-Blend-ref", qBlendRef ) ; // Used for exp
620 
621  fHadronTensorModel = dynamic_cast< const SuSAv2QELHadronTensorModel* >(
622  this->SubAlg("HadronTensorAlg") );
623  assert( fHadronTensorModel );
624 
625  // Load XSec Integrator
626  fXSecIntegrator = dynamic_cast<const XSecIntegratorI *>(
627  this->SubAlg("XSec-Integrator") );
628  assert( fXSecIntegrator );
629 
630  // Fermi momentum tables for scaling
631  this->GetParam( "FermiMomentumTable", fKFTable);
632 
633  // Binding energy lookups for scaling
634  this->GetParam( "RFG-NucRemovalE@Pdg=1000020040", fEbHe );
635  this->GetParam( "RFG-NucRemovalE@Pdg=1000030060", fEbLi );
636  this->GetParam( "RFG-NucRemovalE@Pdg=1000060120", fEbC );
637  this->GetParam( "RFG-NucRemovalE@Pdg=1000080160", fEbO );
638  this->GetParam( "RFG-NucRemovalE@Pdg=1000120240", fEbMg );
639  this->GetParam( "RFG-NucRemovalE@Pdg=1000180400", fEbAr );
640  this->GetParam( "RFG-NucRemovalE@Pdg=1000200400", fEbCa );
641  this->GetParam( "RFG-NucRemovalE@Pdg=1000260560", fEbFe );
642  this->GetParam( "RFG-NucRemovalE@Pdg=1000280580", fEbNi );
643  this->GetParam( "RFG-NucRemovalE@Pdg=1000501190", fEbSn );
644  this->GetParam( "RFG-NucRemovalE@Pdg=1000791970", fEbAu );
645  this->GetParam( "RFG-NucRemovalE@Pdg=1000822080", fEbPb );
646 
647  // Read optional QvalueShifter:
648  fQvalueShifter = nullptr;
649  if( GetConfig().Exists("QvalueShifterAlg") ) {
650 
651  fQvalueShifter = dynamic_cast<const QvalueShifter *> ( this->SubAlg("QvalueShifterAlg") );
652 
653  if( !fQvalueShifter ) {
654 
655  good_config = false ;
656 
657  LOG("SuSAv2QE", pERROR) << "The required QvalueShifterAlg is not valid. AlgID is : "
658  << SubAlg("QvalueShifterAlg")->Id() ;
659  }
660  } // if there is a requested QvalueShifteralgo
661  if( ! good_config ) {
662  LOG("SuSAv2QE", pERROR) << "Configuration has failed.";
663  exit(78) ;
664  }
665 
666 }
#define pERROR
Definition: Messenger.h:59
Cross Section Integrator Interface.
int blendMode
Blending start/end q0 value for the combined models.
const XSecIntegratorI * fXSecIntegrator
GSL numerical integrator.
const QvalueShifter * fQvalueShifter
Creates hadron tensor objects for calculations of quasielastic cross sections using the SuSAv2 approa...
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
double fXSecCCScale
External scaling factor for this cross section.
const HadronTensorModelI * fHadronTensorModel
virtual const AlgId & Id(void) const
Get algorithm ID.
Definition: Algorithm.h:98
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 SuSAv2QELPXSec::ValidProcess ( const Interaction i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 559 of file SuSAv2QELPXSec.cxx.

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

Referenced by XSec().

560 {
561  if ( interaction->TestBit(kISkipProcessChk) ) return true;
562 
563  const InitialState & init_state = interaction->InitState();
564  const ProcessInfo & proc_info = interaction->ProcInfo();
565 
566  if ( !proc_info.IsQuasiElastic() ) return false;
567 
568  // The calculation is only appropriate for complex nuclear targets,
569  // not free nucleons.
570  if ( !init_state.Tgt().IsNucleus() ) return false;
571 
572  int nuc = init_state.Tgt().HitNucPdg();
573  int nu = init_state.ProbePdg();
574 
575  bool isP = pdg::IsProton(nuc);
576  bool isN = pdg::IsNeutron(nuc);
577  bool isnu = pdg::IsNeutrino(nu);
578  bool isnub = pdg::IsAntiNeutrino(nu);
579  bool is_chgl = pdg::IsChargedLepton(nu);
580 
581  bool prcok = ( proc_info.IsWeakCC() && ((isP && isnub) || (isN && isnu)) )
582  || ( proc_info.IsEM() && is_chgl && (isP || isN) );
583  if ( !prcok ) return false;
584 
585  return true;
586 }
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
int HitNucPdg(void) const
Definition: Target.cxx:304
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:69
bool IsNucleus(void) const
Definition: Target.cxx:272
bool IsChargedLepton(int pdgc)
Definition: PDGUtils.cxx:101
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
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:118
int ProbePdg(void) const
Definition: InitialState.h:64
bool IsEM(void) const
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 SuSAv2QELPXSec::XSec ( const Interaction i,
KinePhaseSpace_t  k 
) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 46 of file SuSAv2QELPXSec.cxx.

References genie::KinePhaseSpace::AsString(), blendMode, genie::units::cm2, genie::LabFrameHadronTensorI::dSigma_dT_dCosTheta_rosenbluth(), fEbAr, fEbC, fEbFe, fEbHe, fEbLi, fEbMg, fEbO, fEbPb, fEbSn, fHadronTensorModel, fQvalueShifter, genie::Interaction::FSPrimLepton(), fXSecCCScale, fXSecEMScale, fXSecNCScale, genie::Kinematics::GetKV(), genie::utils::mec::Getq0q3FromTlCostl(), genie::HadronTensorModelI::GetTensor(), genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IonPdgCodeToA(), genie::pdg::IonPdgCodeToZ(), genie::pdg::IsAntiNeutrino(), genie::ProcessInfo::IsEM(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), genie::kHT_QE_CRPA_anu_High, genie::kHT_QE_CRPA_anu_Low, genie::kHT_QE_CRPA_anu_Medium, genie::kHT_QE_CRPA_High, genie::kHT_QE_CRPA_Low, genie::kHT_QE_CRPA_Medium, genie::kHT_QE_CRPAPW_anu_High, genie::kHT_QE_CRPAPW_anu_Low, genie::kHT_QE_CRPAPW_anu_Medium, genie::kHT_QE_CRPAPW_High, genie::kHT_QE_CRPAPW_Low, genie::kHT_QE_CRPAPW_Medium, genie::kHT_QE_EM, genie::kHT_QE_EM_neutron, genie::kHT_QE_EM_proton, genie::kHT_QE_Full, genie::kHT_QE_HF_anu_High, genie::kHT_QE_HF_anu_Low, genie::kHT_QE_HF_anu_Medium, genie::kHT_QE_HF_High, genie::kHT_QE_HF_Low, genie::kHT_QE_HF_Medium, genie::kHT_QE_HFPW_anu_High, genie::kHT_QE_HFPW_anu_Low, genie::kHT_QE_HFPW_anu_Medium, genie::kHT_QE_HFPW_High, genie::kHT_QE_HFPW_Low, genie::kHT_QE_HFPW_Medium, genie::kHT_QE_SuSABlend, genie::kHT_QE_SuSABlend_anu, genie::kHT_Undefined, genie::Interaction::Kine(), genie::kKVctl, genie::kKVTl, kMd_CRPA, kMd_CRPAPW, kMd_CRPAPWSuSAv2Hybrid, kMd_CRPASuSAv2Hybrid, kMd_HF, kMd_HFPW, kMd_HFPWSuSAv2Hybrid, kMd_HFSuSAv2Hybrid, kMd_SuSAv2, kMd_SuSAv2Blend, genie::controls::kMinQ2Limit, genie::kPdgTgtC12, genie::kPdgTgtO16, genie::constants::kPi, genie::kPSTlctl, genie::kRfLab, LOG, modelConfig, pDEBUG, genie::Target::Pdg(), genie::InitialState::ProbeE(), genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), pWARN, q0BlendEnd, q0BlendStart, genie::HadronTensorI::q0Max(), genie::HadronTensorI::q0Min(), genie::utils::kinematics::Q2(), qBlendDel, qBlendRef, genie::HadronTensorI::qMagMax(), genie::HadronTensorI::qMagMin(), genie::utils::mec::Qvalue(), genie::InitialState::Tgt(), ValidProcess(), and XSecScaling().

48 {
49  if ( !this->ValidProcess(interaction) ) return 0.;
50 
51  // Check that the input kinematical point is within the range
52  // in which hadron tensors are known (for chosen target)
53  double Ev = interaction->InitState().ProbeE(kRfLab);
54  double Tl = interaction->Kine().GetKV(kKVTl);
55  double costl = interaction->Kine().GetKV(kKVctl);
56  double ml = interaction->FSPrimLepton()->Mass();
57  double Q0 = 0.;
58  double Q3 = 0.;
59 
60  genie::utils::mec::Getq0q3FromTlCostl(Tl, costl, Ev, ml, Q0, Q3);
61 
62  // *** Enforce the global Q^2 cut (important for EM scattering) ***
63  // Choose the appropriate minimum Q^2 value based on the interaction
64  // mode (this is important for EM interactions since the differential
65  // cross section blows up as Q^2 --> 0)
66  double Q2min = genie::controls::kMinQ2Limit; // CC/NC limit
67  if ( interaction->ProcInfo().IsEM() ) Q2min = genie::utils::kinematics
68  ::electromagnetic::kMinQ2Limit; // EM limit
69 
70  // Neglect shift due to binding energy. The cut is on the actual
71  // value of Q^2, not the effective one to use in the tensor contraction.
72  double Q2 = Q3*Q3 - Q0*Q0;
73  if ( Q2 < Q2min ) return 0.;
74 
75 
76  // ******************************
77  // Now choose which tesor to use
78  // ******************************
79 
80  // Get the hadron tensor for the selected nuclide. Check the probe PDG code
81  // to know whether to use the tensor for CC neutrino scattering or for
82  // electron scattering
83  int target_pdg = interaction->InitState().Tgt().Pdg();
84  int probe_pdg = interaction->InitState().ProbePdg();
85  // get the hit nucleon pdg code
86  int hit_nuc_pdg = interaction->InitState().Tgt().HitNucPdg();
87 
88  int tensor_pdg_susa = target_pdg;
89  int tensor_pdg_crpa = target_pdg;
90  int A_request = pdg::IonPdgCodeToA(target_pdg);
91  int Z_request = pdg::IonPdgCodeToZ(target_pdg);
92  bool need_to_scale_susa = false;
93  bool need_to_scale_crpa = false;
94 
95 
96  // This gets a bit messy as the different models have different
97  // targets available and currently only SuSA does EM
98 
99  HadronTensorType_t tensor_type_susa = kHT_Undefined;
100  HadronTensorType_t tensor_type_crpa = kHT_Undefined;
101  HadronTensorType_t tensor_type_blen = kHT_Undefined;
102 
103  if ( pdg::IsNeutrino(probe_pdg) ) {
104  tensor_type_susa = kHT_QE_Full;
105  tensor_type_blen = kHT_QE_SuSABlend;
106  // CRPA/HF tensors having q0 dependent binning, so are split
107  // CRPA
109  if(Q0<0.060) tensor_type_crpa = kHT_QE_CRPA_Low;
110  else if(Q0<0.150) tensor_type_crpa = kHT_QE_CRPA_Medium;
111  else tensor_type_crpa = kHT_QE_CRPA_High;
112  }
113  // Hartree-Fock
115  if(Q0<0.060) tensor_type_crpa = kHT_QE_HF_Low;
116  else if(Q0<0.150) tensor_type_crpa = kHT_QE_HF_Medium;
117  else tensor_type_crpa = kHT_QE_HF_High;
118  }
120  if(Q0<0.060) tensor_type_crpa = kHT_QE_CRPAPW_Low;
121  else if(Q0<0.150) tensor_type_crpa = kHT_QE_CRPAPW_Medium;
122  else tensor_type_crpa = kHT_QE_CRPAPW_High;
123  }
124  // Hartree-Fock
126  if(Q0<0.060) tensor_type_crpa = kHT_QE_HFPW_Low;
127  else if(Q0<0.150) tensor_type_crpa = kHT_QE_HFPW_Medium;
128  else tensor_type_crpa = kHT_QE_HFPW_High;
129  }
130  }
131  else if ( pdg::IsAntiNeutrino(probe_pdg) ){
132  // SuSA implementation doesn't accoutn for asymmetry between protons
133  // and neutrons. In general this is a small effect.
134  tensor_type_susa = kHT_QE_Full;
135  //For the blending case, Ar40 is treated specially:
136  if(A_request == 40 && Z_request == 18){
137  tensor_type_blen = kHT_QE_SuSABlend_anu;
138  }
139  else tensor_type_blen = kHT_QE_SuSABlend;
140  // CRPA tensors having q0 dependent binning, so are split:
141  //CRPA
143  if(Q0<0.060) tensor_type_crpa = kHT_QE_CRPA_anu_Low;
144  else if(Q0<0.150) tensor_type_crpa = kHT_QE_CRPA_anu_Medium;
145  else tensor_type_crpa = kHT_QE_CRPA_anu_High;
146  }
147  // Hartree-Fock
149  if(Q0<0.060) tensor_type_crpa = kHT_QE_HF_anu_Low;
150  else if(Q0<0.150) tensor_type_crpa = kHT_QE_HF_anu_Medium;
151  else tensor_type_crpa = kHT_QE_HF_anu_High;
152  }
154  if(Q0<0.060) tensor_type_crpa = kHT_QE_CRPAPW_anu_Low;
155  else if(Q0<0.150) tensor_type_crpa = kHT_QE_CRPAPW_anu_Medium;
156  else tensor_type_crpa = kHT_QE_CRPAPW_anu_High;
157  }
158  // Hartree-Fock
160  if(Q0<0.060) tensor_type_crpa = kHT_QE_HFPW_anu_Low;
161  else if(Q0<0.150) tensor_type_crpa = kHT_QE_HFPW_anu_Medium;
162  else tensor_type_crpa = kHT_QE_HFPW_anu_High;
163  }
164  }
165  else {
166  // If the probe is not a neutrino, assume that it's an electron
167  // Currently only avaialble for SuSA. CRPA coming soon(ish)!
168  if ( pdg::IsProton(hit_nuc_pdg) ) {
169  tensor_type_susa = kHT_QE_EM_proton;
170  }
171  else if( pdg::IsNeutron(hit_nuc_pdg) ) {
172  tensor_type_susa = kHT_QE_EM_neutron;
173  }
174  else {
175  tensor_type_susa = kHT_QE_EM; // default
176  }
177 
178  }
179 
180  double Eb_tgt=0;
181  double Eb_ten_susa=0;
182  double Eb_ten_crpa=0;
183 
184  if ( A_request <= 4 ) {
185  // Use carbon tensor for very light nuclei. This is not ideal . . .
186  Eb_tgt = fEbHe;
187  tensor_pdg_susa = kPdgTgtC12;
188  tensor_pdg_crpa = kPdgTgtC12;
189  Eb_ten_susa = fEbC;
190  Eb_ten_crpa = fEbC;
191  }
192  else if (A_request < 9) {
193  Eb_tgt=fEbLi;
194  tensor_pdg_susa = kPdgTgtC12;
195  tensor_pdg_crpa = kPdgTgtC12;
196  Eb_ten_susa = fEbC;
197  Eb_ten_crpa = fEbC;
198  }
199  else if (A_request >= 9 && A_request < 15) {
200  Eb_tgt=fEbC;
201  tensor_pdg_susa = kPdgTgtC12;
202  tensor_pdg_crpa = kPdgTgtC12;
203  Eb_ten_susa = fEbC;
204  Eb_ten_crpa = fEbC;
205  }
206  else if(A_request >= 15 && A_request < 22) {
207  Eb_tgt=fEbO;
208  tensor_pdg_susa = kPdgTgtC12;
209  tensor_pdg_crpa = kPdgTgtO16;
210  Eb_ten_susa = fEbC;
211  Eb_ten_crpa = fEbO;
212  }
213  else if(A_request == 40 && Z_request == 18) {
214  // Treat the common non-isoscalar case specially
215  Eb_tgt=fEbAr;
216  tensor_pdg_susa = kPdgTgtC12;
217  tensor_pdg_crpa = 1000180400;
218  Eb_ten_susa = fEbC;
219  Eb_ten_crpa = fEbAr;
220  }
221  else if(A_request >= 22 && A_request < 40) {
222  Eb_tgt=fEbMg;
223  tensor_pdg_susa = kPdgTgtC12;
224  tensor_pdg_crpa = kPdgTgtO16;
225  Eb_ten_susa = fEbC;
226  Eb_ten_crpa = fEbO;
227  }
228  else if(A_request >= 40 && A_request < 56) {
229  Eb_tgt=fEbAr;
230  tensor_pdg_susa = kPdgTgtC12;
231  tensor_pdg_crpa = kPdgTgtO16;
232  Eb_ten_susa = fEbC;
233  Eb_ten_crpa = fEbO;
234  }
235  else if(A_request >= 56 && A_request < 119) {
236  Eb_tgt=fEbFe;
237  tensor_pdg_susa = kPdgTgtC12;
238  tensor_pdg_crpa = kPdgTgtO16;
239  Eb_ten_susa = fEbC;
240  Eb_ten_crpa = fEbO;
241  }
242  else if(A_request >= 119 && A_request < 206) {
243  Eb_tgt=fEbSn;
244  tensor_pdg_susa = kPdgTgtC12;
245  tensor_pdg_crpa = kPdgTgtO16;
246  Eb_ten_susa = fEbC;
247  Eb_ten_crpa = fEbO;
248  }
249  else if(A_request >= 206) {
250  Eb_tgt=fEbPb;
251  tensor_pdg_susa = kPdgTgtC12;
252  tensor_pdg_crpa = kPdgTgtO16;
253  Eb_ten_susa = fEbC;
254  Eb_ten_crpa = fEbO;
255  }
256 
257  if (tensor_pdg_susa != target_pdg) need_to_scale_susa = true;
258  if (tensor_pdg_crpa != target_pdg) need_to_scale_crpa = true;
259 
260  // Finally we can now get the tensors we need
261 
262  const LabFrameHadronTensorI* tensor_susa;
263  const LabFrameHadronTensorI* tensor_crpa;
264  const LabFrameHadronTensorI* tensor_blen;
265 
266  if( modelConfig == kMd_SuSAv2 ){
267  tensor_susa = dynamic_cast<const LabFrameHadronTensorI*>
268  ( fHadronTensorModel->GetTensor (tensor_pdg_susa, tensor_type_susa) );
269 
270  if ( !tensor_susa ) {
271  LOG("SuSAv2QE", pWARN) << "Failed to load a SuSAv2 hadronic tensor for the"
272  " nuclide " << tensor_pdg_susa;
273  return 0.;
274  }
275  }
276 
280  tensor_blen = dynamic_cast<const LabFrameHadronTensorI*>
281  ( fHadronTensorModel->GetTensor (tensor_pdg_crpa, tensor_type_blen) );
282 
283  // If retrieving the tensor failed, complain and return zero
284  if ( !tensor_blen ) {
285  LOG("SuSAv2QE", pWARN) << "Failed to load a blending SuSAv2 hadronic tensor for the"
286  " nuclide " << tensor_pdg_crpa;
287  return 0.;
288  }
289  }
290 
292 
293  tensor_crpa = dynamic_cast<const LabFrameHadronTensorI*>
294  ( fHadronTensorModel->GetTensor (tensor_pdg_crpa, tensor_type_crpa) );
295 
296  // If retrieving the tensor failed, complain and return zero
297  if ( !tensor_crpa ) {
298  LOG("SuSAv2QE", pWARN) << "Failed to load a CRPA or HF hadronic tensor for the"
299  " nuclide " << tensor_pdg_crpa;
300  return 0.;
301  }
302  }
303 
304  // *****************************
305  // Q_value offset calculation
306  // *****************************
307 
308  // SD: The Q-Value essentially corrects q0 to account for nuclear
309  // binding energy in the Valencia model but this effect is already
310  // in the SuSAv2 and CRPA/HF tensors so I'll set it to 0.
311  // However, if I want to scale I need to account for the altered
312  // binding energy. To first order I can use the Q_value for this
313  double Delta_Q_value_susa = Eb_tgt-Eb_ten_susa;
314  double Delta_Q_value_crpa = Eb_tgt-Eb_ten_crpa;
315  double Delta_Q_value_blen = Eb_tgt-Eb_ten_crpa;
316 
317  // Apply Qvalue relative shift if needed:
318  if( fQvalueShifter ) {
319  // We have the option to add an additional shift on top of the binding energy correction
320  // The QvalueShifter, is a relative shift to the Q_value.
321  // The Q_value was already taken into account in the hadron tensor. Here we recalculate it
322  // to get the right absolute shift.
323  double tensor_Q_value_susa = genie::utils::mec::Qvalue(tensor_pdg_susa,probe_pdg);
324  double total_Q_value_susa = tensor_Q_value_susa + Delta_Q_value_susa ;
325  double Q_value_shift_susa = total_Q_value_susa * fQvalueShifter -> Shift( interaction->InitState().Tgt() ) ;
326 
327  double tensor_Q_value_crpa = genie::utils::mec::Qvalue(tensor_pdg_crpa,probe_pdg);
328  double total_Q_value_crpa = tensor_Q_value_crpa + Delta_Q_value_crpa ;
329  double Q_value_shift_crpa = total_Q_value_crpa * fQvalueShifter -> Shift( interaction->InitState().Tgt() ) ;
330 
331  Delta_Q_value_susa += Q_value_shift_susa;
332  Delta_Q_value_crpa += Q_value_shift_crpa;
333  Delta_Q_value_blen += Q_value_shift_crpa;
334  }
335 
336  // Set the xsec to zero for interactions with q0,q3 outside the requested range
337 
338  if( modelConfig == kMd_SuSAv2){
339  double Q0min = tensor_susa->q0Min();
340  double Q0max = tensor_susa->q0Max();
341  double Q3min = tensor_susa->qMagMin();
342  double Q3max = tensor_susa->qMagMax();
343  if (Q0-Delta_Q_value_susa < Q0min || Q0-Delta_Q_value_susa > Q0max || Q3 < Q3min || Q3 > Q3max) {
344  return 0.0;
345  }
346  }
347 
348  else if ( modelConfig == kMd_CRPA || modelConfig == kMd_HF ||
350  double Q0min = tensor_crpa->q0Min();
351  double Q0max = tensor_crpa->q0Max();
352  double Q3min = tensor_crpa->qMagMin();
353  double Q3max = tensor_crpa->qMagMax();
354  if (Q0-Delta_Q_value_crpa < Q0min || Q0-Delta_Q_value_crpa > Q0max || Q3 < Q3min || Q3 > Q3max) {
355  return 0.0;
356  }
357  }
358 
359  else if ( modelConfig == kMd_SuSAv2Blend){
360  double Q0min = tensor_blen->q0Min();
361  double Q0max = tensor_blen->q0Max();
362  double Q3min = tensor_blen->qMagMin();
363  double Q3max = tensor_blen->qMagMax();
364  if (Q0-Delta_Q_value_blen < Q0min || Q0-Delta_Q_value_blen > Q0max || Q3 < Q3min || Q3 > Q3max) {
365  return 0.0;
366  }
367  }
368 
369  else{ // hybrid (blending) cases. Low kinematics handled by CRPA/HF, high kinematics by blended SuSA
370  double Q0min = tensor_crpa->q0Min();
371  double Q0max = tensor_blen->q0Max();
372  double Q3min = tensor_crpa->qMagMin();
373  double Q3max = tensor_blen->qMagMax();
374  if (Q0-Delta_Q_value_crpa < Q0min || Q0-Delta_Q_value_blen > Q0max || Q3 < Q3min || Q3 > Q3max) {
375  return 0.0;
376  }
377  }
378 
379  // ******************************
380  // Actual xsec calculation
381  // ******************************
382 
383  double xsec_susa = 0;
384  double xsec_crpa = 0;
385  double xsec_blen = 0;
386 
387  if( modelConfig == kMd_SuSAv2 ){
388  // Compute the cross section using the hadron tensor
389  xsec_susa = tensor_susa->dSigma_dT_dCosTheta_rosenbluth(interaction, Delta_Q_value_susa);
390  LOG("SuSAv2QE", pDEBUG) << "SuSAv2 XSec in cm2 / neutron is " << xsec_susa/(units::cm2);
391  xsec_susa = XSecScaling(xsec_susa, interaction, target_pdg, tensor_pdg_susa, need_to_scale_susa);
392  }
393 
394 
398  // Compute the cross section using the hadron tensor
399  xsec_blen = tensor_blen->dSigma_dT_dCosTheta_rosenbluth(interaction, Delta_Q_value_blen);
400  LOG("SuSAv2QE", pDEBUG) << "SuSAv2 (blending) XSec in cm2 / atom is " << xsec_blen/(units::cm2);
401  // The blended SuSAv2 calculation already gives the xsec per atom
402  // For the A-scaling below to make sense we need to transform them to per active nucleon
403  int A_tensor = pdg::IonPdgCodeToA(tensor_pdg_crpa);
404  int Z_tensor = pdg::IonPdgCodeToZ(tensor_pdg_crpa);
405  int N_tensor = A_tensor-Z_tensor;
406 
407  if ( pdg::IsNeutrino(probe_pdg) ) xsec_blen *= 1.0/N_tensor;
408  else if ( pdg::IsAntiNeutrino(probe_pdg) ) xsec_blen *= 1.0/Z_tensor;
409 
410  xsec_blen = XSecScaling(xsec_blen, interaction, target_pdg, tensor_pdg_crpa, need_to_scale_crpa);
411  }
412 
414  // Compute the cross section using the hadron tensor
415  xsec_crpa = tensor_crpa->dSigma_dT_dCosTheta_rosenbluth(interaction, Delta_Q_value_crpa);
416  LOG("SuSAv2QE", pDEBUG) << "CRPA or HF XSec in cm2 / atom is " << xsec_crpa/(units::cm2);
417  // The CRPA calculation already gives the xsec per atom
418  // For the A-scaling below to make sense we need to transform them to per active nucleon
419  int A_tensor = pdg::IonPdgCodeToA(tensor_pdg_crpa);
420  int Z_tensor = pdg::IonPdgCodeToZ(tensor_pdg_crpa);
421  int N_tensor = A_tensor-Z_tensor;
422 
423  if ( pdg::IsNeutrino(probe_pdg) ) xsec_crpa *= 1.0/N_tensor;
424  else if ( pdg::IsAntiNeutrino(probe_pdg) ) xsec_crpa *= 1.0/Z_tensor;
425 
426  // TODO: When we add EM for CRPA need to add how this case it dealt with here
427 
428  xsec_crpa = XSecScaling(xsec_crpa, interaction, target_pdg, tensor_pdg_crpa, need_to_scale_crpa);
429 
430  }
431 
432  // Apply blending if needed
433  double xsec = 0;
434 
435  if( modelConfig == kMd_SuSAv2 ) xsec = xsec_susa;
436  if( modelConfig == kMd_SuSAv2Blend ) xsec = xsec_blen;
437  if( modelConfig == kMd_CRPA || modelConfig == kMd_HF ||
438  modelConfig == kMd_CRPAPW || modelConfig == kMd_HFPW ) xsec = xsec_crpa;
439  else if( modelConfig == kMd_CRPASuSAv2Hybrid ||
442  modelConfig == kMd_HFPWSuSAv2Hybrid ){ // blending cases
443  if(blendMode == 1) // Linear blending in q0
444  if (Q0 < q0BlendStart) xsec = xsec_crpa;
445  else if (Q0 > q0BlendEnd) xsec = xsec_blen;
446  else{
447  double SuSAFrac = (Q0 - q0BlendStart) / (q0BlendEnd - q0BlendStart);
448  double CRPAFrac = 1 - SuSAFrac;
449  xsec = SuSAFrac*xsec_blen + CRPAFrac*xsec_crpa;
450  LOG("SuSAv2QE", pDEBUG) << "Q0 is " << Q0;
451  LOG("SuSAv2QE", pDEBUG) << "SuSAFrac is " << SuSAFrac;
452  LOG("SuSAv2QE", pDEBUG) << "CRPAFrac is " << CRPAFrac;
453  LOG("SuSAv2QE", pDEBUG) << "xsec is " << xsec;
454  }
455  else if(blendMode == 2){ // Exp blending in q (from Alexis)
456  double phi_q = (genie::constants::kPi / 2.) * (1 - 1./(1+std::exp( (Q3 - qBlendRef)/qBlendDel)) );
457  xsec = TMath::Sin(phi_q)*TMath::Sin(phi_q)*xsec_blen + TMath::Cos(phi_q)*TMath::Cos(phi_q)*xsec_crpa;
458  LOG("SuSAv2QE", pDEBUG) << "Q3 is " << Q3;
459  LOG("SuSAv2QE", pDEBUG) << "SuSAFrac is " << TMath::Sin(phi_q)*TMath::Sin(phi_q);
460  LOG("SuSAv2QE", pDEBUG) << "CRPAFrac is " << TMath::Cos(phi_q)*TMath::Cos(phi_q);
461  LOG("SuSAv2QE", pDEBUG) << "xsec is " << xsec;
462  }
463  }
464 
465  // Apply given overall scaling factor
466  double xsec_scale = 1 ;
467  if( interaction->ProcInfo().IsWeakCC() ) xsec_scale = fXSecCCScale;
468  else if( interaction->ProcInfo().IsWeakNC() ) xsec_scale = fXSecNCScale;
469  else if( interaction->ProcInfo().IsEM() ) xsec_scale = fXSecEMScale;
470 
471  xsec *= xsec_scale ;
472 
473  if ( kps != kPSTlctl ) {
474  LOG("SuSAv2QE", pWARN)
475  << "Doesn't support transformation from "
476  << KinePhaseSpace::AsString(kPSTlctl) << " to "
477  << KinePhaseSpace::AsString(kps);
478  xsec = 0.;
479  }
480 
481  return xsec;
482 }
virtual const HadronTensorI * GetTensor(int tensor_pdg, HadronTensorType_t type) const =0
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int blendMode
Blending start/end q0 value for the combined models.
const QvalueShifter * fQvalueShifter
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
virtual double dSigma_dT_dCosTheta_rosenbluth(const Interaction *interaction, double Q_value) const =0
Abstract interface for an object that computes the elements ( , , etc.) and structure functions ( ...
enum genie::HadronTensorType HadronTensorType_t
virtual double q0Max() const =0
static const double kMinQ2Limit
Definition: Controls.h:41
const int kPdgTgtO16
Definition: PDGCodes.h:203
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
double Qvalue(int targetpdg, int nupdg)
Definition: MECUtils.cxx:164
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
double fXSecCCScale
External scaling factor for this cross section.
static constexpr double cm2
Definition: Units.h:69
static string AsString(KinePhaseSpace_t kps)
bool IsAntiNeutrino(int pdgc)
Definition: PDGUtils.cxx:118
bool Getq0q3FromTlCostl(double Tl, double costl, double Enu, double ml, double &q0, double &q3)
Definition: MECUtils.cxx:121
const HadronTensorModelI * fHadronTensorModel
double XSecScaling(double xsec, const Interaction *i, int target_pdg, int tensor_pdg, bool need_to_scale) const
#define pWARN
Definition: Messenger.h:60
virtual double qMagMax() const =0
virtual double q0Min() const =0
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
const int kPdgTgtC12
Definition: PDGCodes.h:202
int IonPdgCodeToZ(int pdgc)
Definition: PDGUtils.cxx:55
virtual double qMagMin() const =0
#define pDEBUG
Definition: Messenger.h:63
double SuSAv2QELPXSec::XSecScaling ( double  xsec,
const Interaction i,
int  target_pdg,
int  tensor_pdg,
bool  need_to_scale 
) const
private

Definition at line 485 of file SuSAv2QELPXSec.cxx.

References genie::ProcessInfo::AsString(), genie::units::cm2, genie::FermiMomentumTable::FindClosestKF(), fKFTable, genie::FermiMomentumTablePool::GetTable(), genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::FermiMomentumTablePool::Instance(), genie::pdg::IsAntiNeutrino(), genie::ProcessInfo::IsEM(), genie::pdg::IsNeutrino(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::ProcessInfo::IsWeakCC(), genie::ProcessInfo::IsWeakNC(), genie::kPdgProton, LOG, genie::Target::N(), pDEBUG, pERROR, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), genie::InitialState::Tgt(), and genie::Target::Z().

Referenced by XSec().

486 {
487  // The xsecs need to be given per active nucleon, but the calculations above are per atom.
488  // We also need to A-scale anyhow if the target nucleus is not exactly the one we have the tensor for.
489  // We adjust for this bellow.
490 
491  const ProcessInfo& proc_info = interaction->ProcInfo();
492 
493  // Neutron, proton, and mass numbers of the target
494  const Target& tgt = interaction->InitState().Tgt();
495 
496  int probe_pdg = interaction->InitState().ProbePdg();
497 
498  if ( proc_info.IsWeakCC() ) {
499  if ( pdg::IsNeutrino(probe_pdg) ) xsec *= tgt.N();
500  else if ( pdg::IsAntiNeutrino(probe_pdg) ) xsec *= tgt.Z();
501  else {
502  // We should never get here if ValidProcess() is working correctly
503  LOG("SuSAv2QE", pERROR) << "Unrecognized probe " << probe_pdg
504  << " encountered for a WeakCC process";
505  xsec = 0.;
506  }
507  }
508  else if ( proc_info.IsEM() || proc_info.IsWeakNC() ) {
509  // For EM processes, scale by the number of nucleons of the same type
510  // as the struck one. This ensures the correct ratio of initial-state
511  // p vs. n when making splines. The nuclear cross section is obtained
512  // by scaling by A/2 for an isoscalar target, so we can get the right
513  // behavior for all targets by scaling by Z/2 or N/2 as appropriate.
514  // Do the same for NC. TODO: double-check that this is the right
515  // thing to do when we SuSAv2 NC hadronic tensors are added to GENIE.
516  int hit_nuc_pdg = tgt.HitNucPdg();
517  if ( pdg::IsProton(hit_nuc_pdg) ) xsec *= tgt.Z() / 2.;
518  else if ( pdg::IsNeutron(hit_nuc_pdg) ) xsec *= tgt.N() / 2.;
519  // We should never get here if ValidProcess() is working correctly
520  else return 0.;
521  }
522  else {
523  // We should never get here if ValidProcess() is working correctly
524  LOG("SuSAv2QE", pERROR) << "Unrecognized process " << proc_info.AsString()
525  << " encountered in SuSAv2QELPXSec::XSec()";
526  xsec = 0.;
527  }
528 
529  LOG("SuSAv2QE", pDEBUG) << "XSec in cm2 / atom is " << xsec / units::cm2;
530 
531  // This scaling should be okay-ish for the total xsec, but it misses
532  // the energy shift. To get this we should really just build releveant
533  // hadron tensors but there may be some ways to approximate it.
534  // For more details see Guille's thesis: https://idus.us.es/xmlui/handle/11441/74826
535 
536  // We already did some of this when we apply the Q-value shift. We can do a little
537  // better by tuning the A-scaling as below, following the SuperScaling ansatz
538  if ( need_to_scale ) {
540  const FermiMomentumTable * kft = kftp->GetTable(fKFTable);
541  double KF_tgt = kft->FindClosestKF(target_pdg, kPdgProton);
542  double KF_ten = kft->FindClosestKF(tensor_pdg, kPdgProton);
543  LOG("SuSAv2QE", pDEBUG) << "KF_tgt = " << KF_tgt;
544  LOG("SuSAv2QE", pDEBUG) << "KF_ten = " << KF_ten;
545  double scaleFact = (KF_ten/KF_tgt); // A-scaling already applied in section above
546  xsec *= scaleFact;
547  }
548 
549  return xsec;
550 }
bool IsWeakCC(void) const
bool IsNeutrino(int pdgc)
Definition: PDGUtils.cxx:110
#define pERROR
Definition: Messenger.h:59
int HitNucPdg(void) const
Definition: Target.cxx:304
static FermiMomentumTablePool * Instance(void)
A table of Fermi momentum constants.
bool IsNeutron(int pdgc)
Definition: PDGUtils.cxx:341
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:336
bool IsWeakNC(void) const
Singleton class to load &amp; 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
const FermiMomentumTable * GetTable(string name)
static constexpr double cm2
Definition: Units.h:69
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:118
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
string AsString(void) const
int Z(void) const
Definition: Target.h:68
bool IsEM(void) const
int N(void) const
Definition: Target.h:69
const int kPdgProton
Definition: PDGCodes.h:81
double FindClosestKF(int target_pdgc, int nucleon_pdgc) const
#define pDEBUG
Definition: Messenger.h:63

Member Data Documentation

int genie::SuSAv2QELPXSec::blendMode
private

Blending start/end q0 value for the combined models.

Definition at line 74 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbAr
private

Definition at line 108 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbAu
private

Definition at line 113 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig().

double genie::SuSAv2QELPXSec::fEbC
private

Definition at line 105 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbCa
private

Definition at line 109 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig().

double genie::SuSAv2QELPXSec::fEbFe
private

Definition at line 110 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbHe
private

Definition at line 103 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbLi
private

Definition at line 104 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbMg
private

Definition at line 107 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbNi
private

Definition at line 111 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig().

double genie::SuSAv2QELPXSec::fEbO
private

Definition at line 106 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbPb
private

Definition at line 114 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fEbSn
private

Definition at line 112 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

const XSecAlgorithmI* genie::SuSAv2QELPXSec::fFreeNucleonXSecAlg
private

Alternate cross section model for free nucleon targets.

Definition at line 120 of file SuSAv2QELPXSec.h.

const HadronTensorModelI* genie::SuSAv2QELPXSec::fHadronTensorModel
private

Definition at line 97 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

string genie::SuSAv2QELPXSec::fKFTable
private

Definition at line 100 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSecScaling().

const QvalueShifter* genie::SuSAv2QELPXSec::fQvalueShifter
private

Definition at line 122 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fXSecCCScale
private

External scaling factor for this cross section.

Definition at line 69 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::fXSecEMScale
private

Definition at line 71 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

const XSecIntegratorI* genie::SuSAv2QELPXSec::fXSecIntegrator
private

GSL numerical integrator.

Definition at line 117 of file SuSAv2QELPXSec.h.

Referenced by Integral(), and LoadConfig().

double genie::SuSAv2QELPXSec::fXSecNCScale
private

Definition at line 70 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

modelType genie::SuSAv2QELPXSec::modelConfig
private

Definition at line 95 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::q0BlendEnd
private

Definition at line 76 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::q0BlendStart
private

Definition at line 75 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::qBlendDel
private

Definition at line 77 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::SuSAv2QELPXSec::qBlendRef
private

Definition at line 78 of file SuSAv2QELPXSec.h.

Referenced by LoadConfig(), and XSec().


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