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

Fit to inelastic cross sections for A(e,e')X valid for all W<3 GeV and all Q2<10 GeV2. More...

#include <BostedChristyEMPXSec.h>

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

Public Member Functions

 BostedChristyEMPXSec ()
 
 BostedChristyEMPXSec (string config)
 
virtual ~BostedChristyEMPXSec ()
 
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...
 
bool ValidKinematics (const Interaction *i) const
 Is the input kinematical point a physically allowed one? More...
 
void Configure (const Registry &config)
 
void Configure (string config)
 
- Public Member Functions inherited from genie::XSecAlgorithmI
virtual ~XSecAlgorithmI ()
 
- Public Member Functions inherited from genie::Algorithm
virtual ~Algorithm ()
 
virtual void FindConfig (void)
 
virtual const RegistryGetConfig (void) const
 
RegistryGetOwnedConfig (void)
 
virtual const AlgIdId (void) const
 Get algorithm ID. More...
 
virtual AlgStatus_t GetStatus (void) const
 Get algorithm status. More...
 
virtual bool AllowReconfig (void) const
 
virtual AlgCmp_t Compare (const Algorithm *alg) const
 Compare with input algorithm. More...
 
virtual void SetId (const AlgId &id)
 Set algorithm ID. More...
 
virtual void SetId (string name, string config)
 
const AlgorithmSubAlg (const RgKey &registry_key) const
 
void AdoptConfig (void)
 
void AdoptSubstructure (void)
 
virtual void Print (ostream &stream) const
 Print algorithm info. More...
 

Private Member Functions

void LoadConfig (void)
 
double sigmaR (int, double, double, bool) const
 
double sigmaNR (int, double, double, bool) const
 
void BranchingRatios (int, double &, double &) const
 
void FermiSmearingD (double, double, double &, double &, double &, double &, bool) const
 
void FermiSmearingA (double, double, double, double, double &, double &, double &, double &) const
 
double FitEMC (double, int) const
 
double MEC2009 (int, double, double) const
 

Private Attributes

bool fUseMEC
 account for MEC contribution? More...
 
double fPM
 mass parameter More...
 
double fMP
 mass parameter More...
 
double fAM
 mass parameter More...
 
double fMD
 deuterium mass More...
 
double fMpi0
 pion mass More...
 
double fMeta
 eta mass More...
 
double fWmin
 minimal W More...
 
double fWmax
 maximal W More...
 
double fQ2min
 minimal Q2 More...
 
double fQ2max
 maximal Q2 More...
 
std::array< std::array< double, 3 >, 7 > fBRp
 branching ratios of resonances for proton fit More...
 
std::array< std::array< double, 3 >, 7 > fBRD
 branching ratios of resonances for deterium fit More...
 
std::array< int, 7 > fAngRes
 resonance angular momentum More...
 
std::array< double, 7 > fMassRes
 resonance mass More...
 
std::array< double, 7 > fWidthRes
 resonance width More...
 
std::array< std::array< double, 4 >, 7 > fRescoefTp
 tunable parameters from Ref.1, Table III for resonance More...
 
std::array< std::array< double, 4 >, 7 > fRescoefTD
 tunable parameters from Ref.2, Table III for resonance More...
 
std::array< std::array< double, 3 >, 7 > fRescoefL
 tunable parameters from Ref.1, Table III for resonance More...
 
std::array< std::array< double, 5 >, 2 > fNRcoefTp
 tunable parameters from Ref.1, Table III for nonres bkg More...
 
std::array< std::array< double, 5 >, 2 > fNRcoefTD
 tunable parameters from Ref.1, Table IV for nonres bkg More...
 
std::array< double, 6 > fNRcoefL
 tunable parameters from Ref.1, Table III for nonres bkg More...
 
std::array< double, 6 > fMECcoef
 tunable parameters for Eqs.(20), (21) Ref.2 More...
 
std::array< double, 8 > fMEC2009coef
 tunable parameters for MEC2009 function More...
 
std::array< double, 13 > fAfitcoef
 tunable parameters for nuclei fit More...
 
std::array< double, 9 > fEMCalpha
 tunable parameters for EMC fit More...
 
std::array< double, 3 > fEMCc
 tunable parameters for EMC fit More...
 
map< int, double > fMEC2009p18
 
map< int, double > fKFTable
 
map< int, double > fNucRmvE
 
const XSecIntegratorIfXSecIntegrator
 

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

Fit to inelastic cross sections for A(e,e')X valid for all W<3 GeV and all Q2<10 GeV2.

Author
Igor Kakorin kakor.nosp@m.in@j.nosp@m.inr.r.nosp@m.u Joint Institute for Nuclear Research based on fortran code provided on Peter Bosted's site: https://userweb.jlab.org/~bosted/fits.html
References:
  1. M.E. Christy, P.E.Bosted, "Empirical fit to precision inclusive electron-proton cross sections in the resonance region", PRC 81 (2010) 055213
    1. P.E.Bosted, M.E.Christy, "Empirical fit to inelastic electron-deuteron and electron-neutron resonance region transverse cross", PRC 77 (2008) 065206
    2. C. Maieron, T. W. Donnelly, and I. Sick, "Extended superscaling of electron scattering from nuclei", PRC 65 (2001) 025502
Created:
April 3, 2021
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 41 of file BostedChristyEMPXSec.h.

Constructor & Destructor Documentation

BostedChristyEMPXSec::BostedChristyEMPXSec ( )

Definition at line 35 of file BostedChristyEMPXSec.cxx.

35  :
36 XSecAlgorithmI("genie::BostedChristyEMPXSec")
37 {
38 }
BostedChristyEMPXSec::BostedChristyEMPXSec ( string  config)

Definition at line 40 of file BostedChristyEMPXSec.cxx.

40  :
41 XSecAlgorithmI("genie::BostedChristyEMPXSec", config)
42 {
43 }
BostedChristyEMPXSec::~BostedChristyEMPXSec ( )
virtual

Definition at line 45 of file BostedChristyEMPXSec.cxx.

46 {
47 
48 }

Member Function Documentation

void BostedChristyEMPXSec::BranchingRatios ( int  respdg,
double &  brpi,
double &  breta 
) const
private

Definition at line 669 of file BostedChristyEMPXSec.cxx.

References genie::PDGLibrary::Find(), genie::PDGLibrary::Instance(), genie::pdg::IsNucleon(), genie::pdg::IsPion(), and genie::kPdgEta.

Referenced by LoadConfig().

670 {
671  brpi = 0.;
672  breta = 0.;
673  PDGLibrary * pdglib = PDGLibrary::Instance();
674  TParticlePDG * res_pdg = pdglib->Find(respdg);
675  if (res_pdg != 0)
676  {
677  for (int nch = 0; nch < res_pdg->NDecayChannels(); nch++)
678  {
679  TDecayChannel * ch = res_pdg->DecayChannel(nch);
680  if (ch->NDaughters() == 2)
681  {
682  int first_daughter_pdg = ch->DaughterPdgCode (0);
683  int second_daughter_pdg = ch->DaughterPdgCode (1);
684  if ((genie::pdg::IsNucleon(first_daughter_pdg ) && genie::pdg::IsPion(second_daughter_pdg)) ||
685  (genie::pdg::IsNucleon(second_daughter_pdg) && genie::pdg::IsPion(first_daughter_pdg )))
686  brpi += ch->BranchingRatio();
687  if (first_daughter_pdg == kPdgEta || second_daughter_pdg == kPdgEta)
688  breta += ch->BranchingRatio();
689  }
690  }
691  }
692 }
bool IsPion(int pdgc)
Definition: PDGUtils.cxx:326
bool IsNucleon(int pdgc)
Definition: PDGUtils.cxx:346
const int kPdgEta
Definition: PDGCodes.h:161
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
void BostedChristyEMPXSec::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 657 of file BostedChristyEMPXSec.cxx.

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

658 {
659  Algorithm::Configure(config);
660  this->LoadConfig();
661 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void BostedChristyEMPXSec::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 663 of file BostedChristyEMPXSec.cxx.

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

664 {
665  Algorithm::Configure(config);
666  this->LoadConfig();
667 }
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
void BostedChristyEMPXSec::FermiSmearingA ( double  Q2,
double  W,
double  pF,
double  Es,
double &  F1p,
double &  F1d,
double &  sigmaT,
double &  sigmaL 
) const
private

Definition at line 199 of file BostedChristyEMPXSec.cxx.

References fMP, fPM, genie::constants::kAem, genie::constants::kPi2, genie::utils::kinematics::Q2(), sigmaNR(), sigmaR(), and genie::utils::kinematics::W().

Referenced by XSec().

200 {
201  // The numbers in arrays bellow were not supposed to change in the original
202  // fortran code and therefore are not configurable
203  static constexpr std::array<double, 99> fyp
204  {0.0272,0.0326,0.0390,0.0464,0.0551,0.0651,0.0766,0.0898,0.1049,
205  0.1221,0.1416,0.1636,0.1883,0.2159,0.2466,0.2807,0.3182,0.3595,
206  0.4045,0.4535,0.5066,0.5637,0.6249,0.6901,0.7593,0.8324,0.9090,
207  0.9890,1.0720,1.1577,1.2454,1.3349,1.4254,1.5163,1.6070,1.6968,
208  1.7849,1.8705,1.9529,2.0313,2.1049,2.1731,2.2350,2.2901,2.3379,
209  2.3776,2.4090,2.4317,2.4454,2.4500,2.4454,2.4317,2.4090,2.3776,
210  2.3379,2.2901,2.2350,2.1731,2.1049,2.0313,1.9529,1.8705,1.7849,
211  1.6968,1.6070,1.5163,1.4254,1.3349,1.2454,1.1577,1.0720,0.9890,
212  0.9090,0.8324,0.7593,0.6901,0.6249,0.5637,0.5066,0.4535,0.4045,
213  0.3595,0.3182,0.2807,0.2466,0.2159,0.1883,0.1636,0.1416,0.1221,
214  0.1049,0.0898,0.0766,0.0651,0.0551,0.0464,0.0390,0.0326,0.0272};
215 
216  static constexpr std::array<double, 99> xxp
217  {-3.000,-2.939,-2.878,-2.816,-2.755,-2.694,-2.633,-2.571,-2.510,
218  -2.449,-2.388,-2.327,-2.265,-2.204,-2.143,-2.082,-2.020,-1.959,
219  -1.898,-1.837,-1.776,-1.714,-1.653,-1.592,-1.531,-1.469,-1.408,
220  -1.347,-1.286,-1.224,-1.163,-1.102,-1.041,-0.980,-0.918,-0.857,
221  -0.796,-0.735,-0.673,-0.612,-0.551,-0.490,-0.429,-0.367,-0.306,
222  -0.245,-0.184,-0.122,-0.061, 0.000, 0.061, 0.122, 0.184, 0.245,
223  0.306, 0.367, 0.429, 0.490, 0.551, 0.612, 0.673, 0.735, 0.796,
224  0.857, 0.918, 0.980, 1.041, 1.102, 1.163, 1.224, 1.286, 1.347,
225  1.408, 1.469, 1.531, 1.592, 1.653, 1.714, 1.776, 1.837, 1.898,
226  1.959, 2.020, 2.082, 2.143, 2.204, 2.265, 2.327, 2.388, 2.449,
227  2.510, 2.571, 2.633, 2.694, 2.755, 2.816, 2.878, 2.939, 3.000};
228 
229  double MN = fPM;
230  double MN2 = MN*MN;
231  double Mp = fMP;
232  double Mp2 = Mp*Mp;
233  double W2 = W*W;
234 
235  double nu = (W2 - MN2 + Q2)/2./MN;
236  double qv = TMath::Sqrt(nu*nu + Q2);
237  // assume this is 2*pf*qv
238  double dW2dpF = 2.*qv;
239  double dW2dEs = 2.*(nu + MN);
240  // switched to using 99 bins!
241  F1p = 0;
242  F1d = 0;
243  sigmaT = 0;
244  sigmaL = 0;
245  for (int i=0; i<99; i++)
246  {
247  double fyuse = fyp[i]/100.;
248  double W2p = W2 + xxp[i]*pF*dW2dpF - Es*dW2dEs;
249  if(W2p>1.159)
250  {
251  //proton
252  double Wp = TMath::Sqrt(W2p);
253  double sigmaTp = sigmaR(0, Q2, Wp) + sigmaNR(0, Q2, Wp);
254  double sigmaLp = sigmaR(1, Q2, Wp) + sigmaNR(1, Q2, Wp);
255  double F1pp = sigmaTp*(W2p-Mp2)/8./kPi2/kAem;
256  //neutron
257  double sigmaTd = sigmaR(0, Q2, Wp, true) + sigmaNR(0, Q2, Wp, true);
258  double F1dp = sigmaTd*(W2p-Mp2)/8./kPi2/kAem;
259  F1d += F1dp*fyuse;
260  F1p += F1pp*fyuse;
261  sigmaT += sigmaTp*fyuse;
262  sigmaL += sigmaLp*fyuse;
263  }
264  }
265 
266 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
double sigmaR(int, double, double, bool) const
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
double sigmaNR(int, double, double, bool) const
void BostedChristyEMPXSec::FermiSmearingD ( double  Q2,
double  W,
double &  F1,
double &  R,
double &  sigmaT,
double &  sigmaL,
bool  isDeuterium = false 
) const
private

Definition at line 269 of file BostedChristyEMPXSec.cxx.

References fAM, fMD, fMECcoef, fMP, fUseMEC, genie::constants::kAem, genie::constants::kPi2, genie::utils::kinematics::Q2(), sigmaNR(), sigmaR(), and genie::utils::kinematics::W().

Referenced by XSec().

270 {
271  // The numbers in arrays bellow were not supposed to change in the original
272  // fortran code and therefore are not configurable
273  static constexpr std::array<double, 20> fyd
274  {0.4965, 0.4988, 0.4958, 0.5008, 0.5027, 0.5041, 0.5029, 0.5034,
275  0.4993, 0.5147, 0.5140, 0.4975, 0.5007, 0.4992, 0.4994, 0.4977,
276  0.5023, 0.4964, 0.4966, 0.4767};
277 
278  static constexpr std::array<double, 20> avpz
279  {-0.1820,-0.0829,-0.0590,-0.0448,-0.0345,-0.0264,-0.0195, -0.0135,
280  -0.0079,-0.0025, 0.0029, 0.0083, 0.0139, 0.0199, 0.0268, 0.0349,
281  0.0453, 0.0598, 0.0844, 0.1853};
282 
283  static constexpr std::array<double, 20> avp2
284  {0.0938, 0.0219, 0.0137, 0.0101, 0.0081, 0.0068, 0.0060, 0.0054,
285  0.0051, 0.0049, 0.0050, 0.0051, 0.0055, 0.0060, 0.0069, 0.0081,
286  0.0102, 0.0140, 0.0225, 0.0964};
287 
288  // Look up tables for deuteron in fine bins for sub threshold
289  static constexpr std::array<double, 200> fydf
290  {0.00001,0.00002,0.00003,0.00005,0.00006,0.00009,0.00010,0.00013,
291  0.00015,0.00019,0.00021,0.00026,0.00029,0.00034,0.00038,0.00044,
292  0.00049,0.00057,0.00062,0.00071,0.00078,0.00089,0.00097,0.00109,
293  0.00119,0.00134,0.00146,0.00161,0.00176,0.00195,0.00211,0.00232,
294  0.00252,0.00276,0.00299,0.00326,0.00352,0.00383,0.00412,0.00447,
295  0.00482,0.00521,0.00560,0.00603,0.00648,0.00698,0.00747,0.00803,
296  0.00859,0.00921,0.00985,0.01056,0.01126,0.01205,0.01286,0.01376,
297  0.01467,0.01569,0.01671,0.01793,0.01912,0.02049,0.02196,0.02356,
298  0.02525,0.02723,0.02939,0.03179,0.03453,0.03764,0.04116,0.04533,
299  0.05004,0.05565,0.06232,0.07015,0.07965,0.09093,0.10486,0.12185,
300  0.14268,0.16860,0.20074,0.24129,0.29201,0.35713,0.44012,0.54757,
301  0.68665,0.86965,1.11199,1.43242,1.86532,2.44703,3.22681,4.24972,
302  5.54382,7.04016,8.48123,9.40627,9.40627,8.48123,7.04016,5.54382,
303  4.24972,3.22681,2.44703,1.86532,1.43242,1.11199,0.86965,0.68665,
304  0.54757,0.44012,0.35713,0.29201,0.24129,0.20074,0.16860,0.14268,
305  0.12185,0.10486,0.09093,0.07965,0.07015,0.06232,0.05565,0.05004,
306  0.04533,0.04116,0.03764,0.03453,0.03179,0.02939,0.02723,0.02525,
307  0.02356,0.02196,0.02049,0.01912,0.01793,0.01671,0.01569,0.01467,
308  0.01376,0.01286,0.01205,0.01126,0.01056,0.00985,0.00921,0.00859,
309  0.00803,0.00747,0.00698,0.00648,0.00603,0.00560,0.00521,0.00482,
310  0.00447,0.00412,0.00383,0.00352,0.00326,0.00299,0.00276,0.00252,
311  0.00232,0.00211,0.00195,0.00176,0.00161,0.00146,0.00134,0.00119,
312  0.00109,0.00097,0.00089,0.00078,0.00071,0.00062,0.00057,0.00049,
313  0.00044,0.00038,0.00034,0.00029,0.00026,0.00021,0.00019,0.00015,
314  0.00013,0.00010,0.00009,0.00006,0.00005,0.00003,0.00002,0.00001};
315 
316  static constexpr std::array<double, 200> avp2f
317  {1.0,0.98974,0.96975,0.96768,0.94782,0.94450,0.92494,0.92047,
318  0.90090,0.89563,0.87644,0.87018,0.85145,0.84434,0.82593,0.81841,
319  0.80021,0.79212,0.77444,0.76553,0.74866,0.73945,0.72264,0.71343,
320  0.69703,0.68740,0.67149,0.66182,0.64631,0.63630,0.62125,0.61154,
321  0.59671,0.58686,0.57241,0.56283,0.54866,0.53889,0.52528,0.51581,
322  0.50236,0.49291,0.47997,0.47063,0.45803,0.44867,0.43665,0.42744,
323  0.41554,0.40656,0.39511,0.38589,0.37488,0.36611,0.35516,0.34647,
324  0.33571,0.32704,0.31656,0.30783,0.29741,0.28870,0.27820,0.26945,
325  0.25898,0.25010,0.23945,0.23023,0.21943,0.20999,0.19891,0.18911,
326  0.17795,0.16793,0.15669,0.14667,0.13553,0.12569,0.11504,0.10550,
327  0.09557,0.08674,0.07774,0.06974,0.06184,0.05484,0.04802,0.04203,
328  0.03629,0.03129,0.02654,0.02247,0.01867,0.01545,0.01251,0.01015,
329  0.00810,0.00664,0.00541,0.00512,0.00512,0.00541,0.00664,0.00810,
330  0.01015,0.01251,0.01545,0.01867,0.02247,0.02654,0.03129,0.03629,
331  0.04203,0.04802,0.05484,0.06184,0.06974,0.07774,0.08674,0.09557,
332  0.10550,0.11504,0.12569,0.13553,0.14667,0.15669,0.16793,0.17795,
333  0.18911,0.19891,0.20999,0.21943,0.23023,0.23945,0.25010,0.25898,
334  0.26945,0.27820,0.28870,0.29741,0.30783,0.31656,0.32704,0.33571,
335  0.34647,0.35516,0.36611,0.37488,0.38589,0.39511,0.40656,0.41554,
336  0.42744,0.43665,0.44867,0.45803,0.47063,0.47997,0.49291,0.50236,
337  0.51581,0.52528,0.53889,0.54866,0.56283,0.57241,0.58686,0.59671,
338  0.61154,0.62125,0.63630,0.64631,0.66182,0.67149,0.68740,0.69703,
339  0.71343,0.72264,0.73945,0.74866,0.76553,0.77444,0.79212,0.80021,
340  0.81841,0.82593,0.84434,0.85145,0.87018,0.87644,0.89563,0.90090,
341  0.92047,0.92494,0.94450,0.94782,0.96768,0.96975,0.98974,1.0};
342 
343  double W2=W*W;
344  double MN = fAM;
345  double MN2 = MN*MN;
346  double MD = fMD;
347  double Mp = fMP;
348  double Mp2 = Mp*Mp;
349  double nu = (W2 - MN2 + Q2)/2./MN;
350  double qv = TMath::Sqrt(nu*nu + Q2);
351  F1 = 0.;
352  R = 0.;
353  sigmaT = 0.;
354  sigmaL = 0.;
355  // Do fast 20 bins if abvoe threshold
356  if(W2>1.30)
357  {
358  for (int ism = 0; ism<20; ism++)
359  {
360  double W2p = TMath::Power(MD + nu - sqrt(MN2 + avp2[ism]),2) - qv*qv + 2.*qv*avpz[ism] - avp2[ism];
361  if(W2p>1.155)
362  {
363  double Wp = TMath::Sqrt(W2p);
364  double sigtp = sigmaR(0, Q2, Wp, isDeuterium) + sigmaNR(0, Q2, Wp, isDeuterium);
365  double F1p = sigtp*(W2p-Mp2)/8./kPi2/kAem;
366  F1 += F1p*fyd[ism]/10.;
367  if (!isDeuterium)
368  {
369  double siglp = sigmaR(1, Q2, Wp) + sigmaNR(1, Q2, Wp);
370  sigmaL += siglp*fyd[ism]/10.;
371  sigmaT += sigtp*fyd[ism]/10.;
372  }
373  }
374  }
375  }
376  else
377  {
378  for (int ism = 0;ism<200;ism++)
379  {
380  double pz = -1. + 0.01*(ism + 0.5);
381  // Need avp2f term to get right behavior x>1!
382  double W2p = TMath::Power(MD + nu - sqrt(MN2 + avp2f[ism]),2) - qv*qv + 2.*qv*pz - avp2f[ism];
383  if(W2p>1.155)
384  {
385  double Wp = TMath::Sqrt(W2p);
386  double sigtp = sigmaR(0, Q2, Wp, isDeuterium) + sigmaNR(0, Q2, Wp, isDeuterium);
387  double F1p = sigtp*(W2p-Mp2)/8./kPi2/kAem;
388  F1 += F1p*fydf[ism]/100.;
389  if (!isDeuterium)
390  {
391  double siglp = sigmaR(1, Q2, Wp) + sigmaNR(1, Q2, Wp);
392  sigmaT += sigtp*fydf[ism]/100.;
393  sigmaL += siglp*fydf[ism]/100.;
394  }
395  }
396  }
397  }
398  if (isDeuterium && fUseMEC)
399  // Ref.2, Eq. (20)
400  F1 += fMECcoef[0]*TMath::Exp(-(W - fMECcoef[1])*(W - fMECcoef[1])/fMECcoef[2])/
401  TMath::Power(1. + TMath::Max(0.3,Q2)/fMECcoef[3],fMECcoef[4])*TMath::Power(nu, fMECcoef[5]);
402  if(!isDeuterium && sigmaT!=0.)
403  R = sigmaL/sigmaT;
404 
405 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
double sigmaR(int, double, double, bool) const
std::array< double, 6 > fMECcoef
tunable parameters for Eqs.(20), (21) Ref.2
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
bool fUseMEC
account for MEC contribution?
double sigmaNR(int, double, double, bool) const
double BostedChristyEMPXSec::FitEMC ( double  x,
int  A 
) const
private

Definition at line 582 of file BostedChristyEMPXSec.cxx.

References fEMCalpha, and fEMCc.

Referenced by XSec().

583 {
584  double fitemc = 1.;
585  if(A<=2)
586  return fitemc;
587 
588  double x_u;
589  if (x>0.70 || x<0.0085)
590  //Out of range of fit
591  {
592  if(x<0.0085)
593  x_u = .0085;
594  if(x>0.70)
595  x_u = 0.70;
596  }
597  else
598  x_u = x;
599 
600  double ln_c = fEMCc[0];
601  for (int i = 1; i<=2; i++)
602  ln_c += fEMCc[i]*TMath::Power(TMath::Log(x_u), i);
603  double c = TMath::Exp(ln_c);
604 
605  double alpha = fEMCalpha[0];
606  for (int i = 1; i<=8; i++)
607  alpha += fEMCalpha[i]*TMath::Power(x_u, i);
608 
609  fitemc = c*TMath::Power(A, alpha);
610  return fitemc;
611 }
std::array< double, 9 > fEMCalpha
tunable parameters for EMC fit
std::array< double, 3 > fEMCc
tunable parameters for EMC fit
static constexpr double A
Definition: Units.h:74
double BostedChristyEMPXSec::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 613 of file BostedChristyEMPXSec.cxx.

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

614 {
615  double xsec = fXSecIntegrator->Integrate(this,interaction);
616  return xsec;
617 }
const XSecIntegratorI * fXSecIntegrator
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
void BostedChristyEMPXSec::LoadConfig ( void  )
private

Definition at line 694 of file BostedChristyEMPXSec.cxx.

References genie::units::A, BranchingRatios(), fAfitcoef, fAM, fAngRes, fBRD, fBRp, fEMCalpha, fEMCc, genie::PDGLibrary::Find(), fKFTable, fMassRes, fMD, fMEC2009coef, fMEC2009p18, fMECcoef, fMeta, fMP, fMpi0, fNRcoefL, fNRcoefTD, fNRcoefTp, fNucRmvE, fPM, fQ2max, fQ2min, fRescoefL, fRescoefTD, fRescoefTp, genie::utils::res::FromPdgCode(), fUseMEC, fWidthRes, fWmax, fWmin, genie::Algorithm::GetConfig(), genie::Registry::GetItemMap(), genie::Algorithm::GetParam(), genie::Algorithm::GetParamDef(), genie::Algorithm::GetParamVect(), genie::PDGLibrary::Instance(), genie::pdg::IonPdgCodeToA(), genie::kPdgD13m1520_N0, genie::kPdgD13m1520_NP, genie::kPdgEta, genie::kPdgF15m1680_N0, genie::kPdgF15m1680_NP, genie::kPdgF37m1950_Delta0, genie::kPdgF37m1950_DeltaM, genie::kPdgF37m1950_DeltaP, genie::kPdgF37m1950_DeltaPP, genie::kPdgP11m1440_N0, genie::kPdgP11m1440_NP, genie::kPdgP33m1232_Delta0, genie::kPdgP33m1232_DeltaM, genie::kPdgP33m1232_DeltaP, genie::kPdgP33m1232_DeltaPP, genie::kPdgPi0, genie::kPdgProton, genie::kPdgS11m1535_N0, genie::kPdgS11m1535_NP, genie::kPdgS11m1650_N0, genie::kPdgS11m1650_NP, genie::kPdgTgtDeuterium, LOG, genie::utils::res::Mass(), genie::utils::res::OrbitalAngularMom(), pALERT, pFATAL, and genie::utils::res::Width().

Referenced by Configure().

695 {
696 
697  PDGLibrary * pdglib = PDGLibrary::Instance();
698  GetParamDef("BostedChristyFitEM-PM", fPM, pdglib->Find(kPdgProton)->Mass());
699  GetParamDef("BostedChristyFitEM-MP", fMP, pdglib->Find(kPdgProton)->Mass());
700  GetParamDef("BostedChristyFitEM-AM", fAM, pdglib->Find(kPdgProton)->Mass());
701  GetParamDef("BostedChristyFitEM-MD", fMD, pdglib->Find(kPdgTgtDeuterium)->Mass());
702  GetParamDef("BostedChristyFitEM-Mpi0", fMpi0, pdglib->Find(kPdgPi0)->Mass());
703  GetParamDef("BostedChristyFitEM-Meta", fMeta, pdglib->Find(kPdgEta)->Mass());
704  GetParamDef("BostedChristyFitEM-Wmin", fWmin, 0.0);
705  GetParamDef("BostedChristyFitEM-Wmax", fWmax, 3.0);
706  GetParamDef("BostedChristyFitEM-Q2min", fQ2min, 0.0);
707  GetParamDef("BostedChristyFitEM-Q2max", fQ2max, 10.0);
708  GetParamDef("BostedChristyFitEM-UseMEC", fUseMEC, true);
709 
710  double BRpi, BReta;
711  double brpi, breta;
712 
713  std::vector<double> vBRpi1;
714  std::vector<double> vBRpi2;
715  std::vector<double> vBReta;
716 
717  // load braching ratios for pi
718  bool useBRpi1Default = (GetParamVect("BostedChristyFitEM-PionBRp", vBRpi1, false)<7);
719  bool useBRpi2Default = (GetParamVect("BostedChristyFitEM-PionBRD", vBRpi2, false)<7);
720  // load braching ratios for eta
721  bool useBRetaDefault = (GetParamVect("BostedChristyFitEM-EtaBR", vBReta, false)<7);
722 
723  if (useBRpi1Default || useBRpi2Default || useBRetaDefault)
724  {
725  // use default branching ratios from PDG table
726  // P33(1232)
727  BRpi = 0.;
728  BReta = 0.;
729  BranchingRatios(kPdgP33m1232_DeltaM, brpi, breta);
730  BRpi += brpi;
731  BReta += breta;
732  BranchingRatios(kPdgP33m1232_Delta0, brpi, breta);
733  BRpi += brpi;
734  BReta += breta;
735  BranchingRatios(kPdgP33m1232_DeltaP, brpi, breta);
736  BRpi += brpi;
737  BReta += breta;
739  BRpi += brpi;
740  BReta += breta;
741  BRpi /= 4.;
742  BReta /= 4.;
743  fBRp[0][0] = BRpi;
744  fBRp[0][2] = BReta;
745  fBRD[0][0] = BRpi;
746 
747  // S11(1535)
748  BRpi = 0.;
749  BReta = 0.;
750  BranchingRatios(kPdgS11m1535_N0, brpi, breta);
751  BRpi += brpi;
752  BReta += breta;
753  BranchingRatios(kPdgS11m1535_NP, brpi, breta);
754  BRpi += brpi;
755  BReta += breta;
756  BRpi /= 2.;
757  BReta /= 2.;
758  fBRp[1][0] = BRpi;
759  fBRp[1][2] = BReta;
760  fBRD[1][0] = BRpi;
761 
762  // D13(1520)
763  BRpi = 0.;
764  BReta = 0.;
765  BranchingRatios(kPdgD13m1520_N0, brpi, breta);
766  BRpi += brpi;
767  BReta += breta;
768  BranchingRatios(kPdgD13m1520_NP, brpi, breta);
769  BRpi += brpi;
770  BReta += breta;
771  BRpi /= 2.;
772  BReta /= 2.;
773  fBRp[2][0] = BRpi;
774  fBRp[2][2] = BReta;
775  fBRD[2][0] = BRpi;
776 
777  // F15(1680)
778  BRpi = 0.;
779  BReta = 0.;
780  BranchingRatios(kPdgF15m1680_N0, brpi, breta);
781  BRpi += brpi;
782  BReta += breta;
783  BranchingRatios(kPdgF15m1680_NP, brpi, breta);
784  BRpi += brpi;
785  BReta += breta;
786  BRpi /= 2.;
787  BReta /= 2.;
788  fBRp[3][0] = BRpi;
789  fBRp[3][2] = BReta;
790  fBRD[3][0] = BRpi;
791 
792  // S11(1650)
793  BRpi = 0.;
794  BReta = 0.;
795  BranchingRatios(kPdgS11m1650_N0, brpi, breta);
796  BRpi += brpi;
797  BReta += breta;
798  BranchingRatios(kPdgS11m1650_NP, brpi, breta);
799  BRpi += brpi;
800  BReta += breta;
801  BRpi /= 2.;
802  BReta /= 2.;
803  fBRp[4][0] = BRpi;
804  fBRp[4][2] = BReta;
805  fBRD[4][0] = BRpi;
806 
807  // P11(1440)
808  BRpi = 0.;
809  BReta = 0.;
810  BranchingRatios(kPdgP11m1440_N0, brpi, breta);
811  BRpi += brpi;
812  BReta += breta;
813  BranchingRatios(kPdgP11m1440_NP, brpi, breta);
814  BRpi += brpi;
815  BReta += breta;
816  BRpi /= 2.;
817  BReta /= 2.;
818  fBRp[5][0] = BRpi;
819  fBRp[5][2] = BReta;
820  fBRD[5][0] = BRpi;
821 
822  // F37(1950)
823  BRpi = 0.;
824  BReta = 0.;
825  BranchingRatios(kPdgF37m1950_DeltaM , brpi, breta);
826  BRpi += brpi;
827  BReta += breta;
828  BranchingRatios(kPdgF37m1950_Delta0, brpi, breta);
829  BRpi += brpi;
830  BReta += breta;
831  BranchingRatios(kPdgF37m1950_DeltaP, brpi, breta);
832  BRpi += brpi;
833  BReta += breta;
835  BRpi += brpi;
836  BReta += breta;
837  BRpi /= 4.;
838  BReta /= 4.;
839  fBRp[6][0] = BRpi;
840  fBRp[6][2] = BReta;
841  fBRD[6][0] = BRpi;
842  }
843  if (!useBRpi1Default)
844  // single pion branching ratios from config file
845  for (int i=0; i<7; i++)
846  fBRp[i][0] = vBRpi1[i];
847  if (!useBRpi2Default)
848  // single pion branching ratios from config file
849  for (int i=0; i<7; i++)
850  fBRD[i][0] = vBRpi2[i];
851  if (!useBRetaDefault)
852  // eta branching ratios from config file
853  for (int i=0; i<7; i++)
854  fBRp[i][2] = vBReta[i];
855 
856  for (int i=0; i<7; i++)
857  fBRD[i][2] = 0.;
858 
859  if (useBRpi1Default || useBRpi2Default)
860  LOG("BostedChristyEMPXSec", pALERT) << "*** Use branching ratios for pion decay from PDG table";
861 
862  if (useBRetaDefault)
863  LOG("BostedChristyEMPXSec", pALERT) << "*** Use branching ratios for eta decay from PDG table";
864 
865  // 2-pion branching ratios
866  for (int i=0;i<7;i++)
867  {
868  fBRp[i][1] = 1.-fBRp[i][0]-fBRp[i][2];
869  fBRD[i][1] = 1.-fBRD[i][0]-fBRD[i][2];
870  }
871 
872  // Meson angular momentum
880 
881  std::vector<double> vResMass;
882  // load resonance masses
883  bool useResMassDefault = (GetParamVect("BostedChristyFitEM-ResMass", vResMass, false)<7);
884 
885  if (useResMassDefault)
886  {
887  LOG("BostedChristyEMPXSec", pALERT) << "*** Use resonance masses from PDG table";
888  // Resonance mass
894  fMassRes[5] = res::Mass(res::FromPdgCode(kPdgP11m1440_N0)); // P11(1440) roper
896  }
897  else
898  // eta branching ratios from config file
899  for (int i=0; i<7; i++)
900  fMassRes[i] = vResMass[i];
901 
902  std::vector<double> vResWidth;
903  // load resonance masses
904  bool useResWidthDefault = (GetParamVect("BostedChristyFitEM-ResWidth", vResWidth, false)<7);
905 
906  if (useResWidthDefault)
907  {
908  LOG("BostedChristyEMPXSec", pALERT) << "*** Use resonance widths from PDG table";
909  // Resonance width
915  fWidthRes[5] = res::Width(res::FromPdgCode(kPdgP11m1440_N0)); // P11(1440) roper
917  }
918  else
919  // eta branching ratios from config file
920  for (int i=0; i<7; i++)
921  fWidthRes[i] = vResWidth[i];
922 
923  int length;
924 
925  std::vector<double> vRescoef;
926  length = 7;
927  bool isOk = (GetParamVect("BostedChristyFitEM-ResAT0p", vRescoef)>=length);
928  if (!isOk)
929  {
930  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton AT(0)-parameters for xsec^R_T in the config file!";
931  exit(1);
932  }
933  // Ref.1, Table III, AT(0)
934  for (int i=0;i<length;i++)
935  fRescoefTp[i][0] = vRescoef[i];
936 
937  length = 7;
938  isOk = (GetParamVect("BostedChristyFitEM-Resap", vRescoef)>=length);
939  if (!isOk)
940  {
941  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton a-parameters for xsec^R_T in the config file!";
942  exit(1);
943  }
944  // Ref.1, Table III, a
945  for (int i=0;i<length;i++)
946  fRescoefTp[i][1] = vRescoef[i];
947 
948  length = 7;
949  isOk = (GetParamVect("BostedChristyFitEM-Resbp", vRescoef)>=length);
950  if (!isOk)
951  {
952  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton b-parameters parameters for xsec^R_T in the config file!";
953  exit(1);
954  }
955  // Ref.1, Table III, b
956  for (int i=0;i<length;i++)
957  fRescoefTp[i][2] = vRescoef[i];
958 
959  length = 7;
960  isOk = (GetParamVect("BostedChristyFitEM-Rescp", vRescoef)>=length);
961  if (!isOk)
962  {
963  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton c-parameters parameters for xsec^R_T in the config file!";
964  exit(1);
965  }
966  // Ref.1, Table III, c
967  for (int i=0;i<length;i++)
968  fRescoefTp[i][3] = vRescoef[i];
969 
970  length = 7;
971  isOk = (GetParamVect("BostedChristyFitEM-ResAT0D", vRescoef)>=length);
972  if (!isOk)
973  {
974  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium AT(0)-parameters for xsec^R_T in the config file!";
975  exit(1);
976  }
977  // Ref.2, Table III, AT(0)
978  for (int i=0;i<length;i++)
979  fRescoefTD[i][0] = vRescoef[i];
980 
981  length = 7;
982  isOk = (GetParamVect("BostedChristyFitEM-ResaD", vRescoef)>=length);
983  if (!isOk)
984  {
985  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium a-parameters for xsec^R_T in the config file!";
986  exit(1);
987  }
988  // Ref.2, Table III, a
989  for (int i=0;i<length;i++)
990  fRescoefTD[i][1] = vRescoef[i];
991 
992  length = 7;
993  isOk = (GetParamVect("BostedChristyFitEM-ResbD", vRescoef)>=length);
994  if (!isOk)
995  {
996  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium b-parameters parameters for xsec^R_T in the config file!";
997  exit(1);
998  }
999  // Ref.2, Table III, b
1000  for (int i=0;i<length;i++)
1001  fRescoefTD[i][2] = vRescoef[i];
1002 
1003  length = 7;
1004  isOk = (GetParamVect("BostedChristyFitEM-RescD", vRescoef)>=length);
1005  if (!isOk)
1006  {
1007  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium c-parameters parameters for xsec^R_T in the config file!";
1008  exit(1);
1009  }
1010  // Ref.2, Table III, c
1011  for (int i=0;i<length;i++)
1012  fRescoefTD[i][3] = vRescoef[i];
1013 
1014  length = 7;
1015  isOk = (GetParamVect("BostedChristyFitEM-ResAL0", vRescoef)>=length);
1016  if (!isOk)
1017  {
1018  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton AL0-parameters parameters for xsec^R_T in the config file!";
1019  exit(1);
1020  }
1021  // Ref.1, Table III, AL(0)
1022  for (int i=0;i<length;i++)
1023  fRescoefL[i][0] = vRescoef[i];
1024 
1025  length = 7;
1026  isOk = (GetParamVect("BostedChristyFitEM-Resd", vRescoef)>=length);
1027  if (!isOk)
1028  {
1029  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton d-parameters parameters for xsec^R_L in the config file!";
1030  exit(1);
1031  }
1032  // Ref.1, Table III, d
1033  for (int i=0;i<length;i++)
1034  fRescoefL[i][1] = vRescoef[i];
1035 
1036  length = 7;
1037  isOk = (GetParamVect("BostedChristyFitEM-Rese", vRescoef)>=length);
1038  if (!isOk)
1039  {
1040  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton e-parameters parameters for xsec^R_L in the config file!";
1041  exit(1);
1042  }
1043  // Ref.1, Table III, e
1044  for (int i=0;i<length;i++)
1045  fRescoefL[i][2] = vRescoef[i];
1046 
1047 
1048  std::vector<double> vNRcoef;
1049  length = 5;
1050  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT1p", vNRcoef)>=length);
1051  if (!isOk)
1052  {
1053  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton bkg parameters for xsec^NR_T in the config file!";
1054  exit(1);
1055  }
1056  // Ref.1, Table IV: \sigma^NR,1_T(0), aT_1, bT_1, cT_1, dT_1
1057  for (int i=0;i<length;i++)
1058  fNRcoefTp[0][i] = vNRcoef[i];
1059 
1060  length = 5;
1061  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT2p", vNRcoef)>=length);
1062  if (!isOk)
1063  {
1064  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton bkg parameters for xsec^NR_T in the config file!";
1065  exit(1);
1066  }
1067  // Ref.1, Table IV: \sigma^NR,2_T(0), aT_2, bT_2, cT_2, dT_2
1068  for (int i=0;i<length;i++)
1069  fNRcoefTp[1][i] = vNRcoef[i];
1070 
1071  length = 5;
1072  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT1D", vNRcoef)>=length);
1073  if (!isOk)
1074  {
1075  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium bkg parameters for xsec^NR_T in the config file!";
1076  exit(1);
1077  }
1078  // Ref.2, Table IV: \sigma^NR,1_T(0), aT_1, bT_1, cT_1, dT_1
1079  for (int i=0;i<length;i++)
1080  fNRcoefTD[0][i] = vNRcoef[i];
1081 
1082  length = 5;
1083  isOk = (GetParamVect("BostedChristyFitEM-NRXSecT2D", vNRcoef)>=length);
1084  if (!isOk)
1085  {
1086  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough deuterium bkg parameters for xsec^NR_T in the config file!";
1087  exit(1);
1088  }
1089  // Ref.2, Table IV: \sigma^NR,2_T(0), aT_2, bT_2, cT_2, dT_2
1090  for (int i=0;i<length;i++)
1091  fNRcoefTD[1][i] = vNRcoef[i];
1092 
1093  length = 6;
1094  isOk = (GetParamVect("BostedChristyFitEM-NRXSecL", vNRcoef)>=length);
1095  if (!isOk)
1096  {
1097  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough proton bkg parameters for xsec^NR_L in the config file!";
1098  exit(1);
1099  }
1100  // Ref.1, Table IV: \sigma^NR_L, aL, bL, cL, dL, eL
1101  for (int i=0;i<length;i++)
1102  fNRcoefL[i] = vNRcoef[i];
1103 
1104  length = 6;
1105  isOk = (GetParamVect("BostedChristyFitEM-MEC", vNRcoef)>=length);
1106  if (!isOk)
1107  {
1108  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough parameters for MEC in the config file!";
1109  exit(1);
1110  }
1111  for (int i=0;i<length;i++)
1112  fMECcoef[i] = vNRcoef[i];
1113 
1114  length = 8;
1115  isOk = (GetParamVect("BostedChristyFitEM-MEC2009", vNRcoef)>=length);
1116  if (!isOk)
1117  {
1118  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough parameters for MEC2009 in the config file!";
1119  exit(1);
1120  }
1121  for (int i=0;i<length;i++)
1122  fMEC2009coef[i] = vNRcoef[i];
1123 
1124  length = 13;
1125  isOk = (GetParamVect("BostedChristyFitEM-Afit", vNRcoef)>=length);
1126  if (!isOk)
1127  {
1128  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough parameters for nuclei fit (A-fit) in the config file!";
1129  exit(1);
1130  }
1131  for (int i=0;i<length;i++)
1132  fAfitcoef[i] = vNRcoef[i];
1133 
1134 
1135  length = 9;
1136  isOk = (GetParamVect("BostedChristyFitEM-EMCalpha", vNRcoef)>=length);
1137  if (!isOk)
1138  {
1139  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough alpha coefficients for EMC correction in the config file!";
1140  exit(1);
1141  }
1142  for (int i=0;i<length;i++)
1143  fEMCalpha[i] = vNRcoef[i];
1144 
1145  length = 3;
1146  isOk = (GetParamVect("BostedChristyFitEM-EMCc", vNRcoef)>=length);
1147  if (!isOk)
1148  {
1149  LOG("BostedChristyEMPXSec", pFATAL) << "*** Can't find enough c coefficients for EMC correction in the config file!";
1150  exit(1);
1151  }
1152  for (int i=0;i<length;i++)
1153  fEMCc[i] = vNRcoef[i];
1154 
1155 
1156  std::string keyStart = "BostedChristy-SeparationE@Pdg=";
1157  RgIMap entries = GetConfig().GetItemMap();
1158  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it)
1159  {
1160  const std::string& key = it->first;
1161  int pdg = 0;
1162  int A = 0;
1163  if (0 == key.compare(0, keyStart.size(), keyStart.c_str()))
1164  {
1165  pdg = atoi(key.c_str() + keyStart.size());
1166  A = pdg::IonPdgCodeToA(pdg);
1167  }
1168  if (0 != pdg && A != 0)
1169  {
1170  std::ostringstream key_ss ;
1171  key_ss << keyStart << pdg;
1172  RgKey rgkey = key_ss.str();
1173  double eb;
1174  GetParam( rgkey, eb) ;
1175  eb = TMath::Max(eb, 0.);
1176  fNucRmvE.insert(map<int,double>::value_type(A,eb));
1177  }
1178  }
1179 
1180  keyStart = "BostedChristy-FermiMomentum@Pdg=";
1181  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it)
1182  {
1183  const std::string& key = it->first;
1184  int pdg = 0;
1185  int A = 0;
1186  if (0 == key.compare(0, keyStart.size(), keyStart.c_str()))
1187  {
1188  pdg = atoi(key.c_str() + keyStart.size());
1189  A = pdg::IonPdgCodeToA(pdg);
1190  }
1191  if (0 != pdg && A != 0)
1192  {
1193  std::ostringstream key_ss ;
1194  key_ss << keyStart << pdg;
1195  RgKey rgkey = key_ss.str();
1196  double pf;
1197  GetParam( rgkey, pf) ;
1198  pf = TMath::Max(pf, 0.);
1199  fKFTable.insert(map<int,double>::value_type(A,pf));
1200  }
1201  }
1202 
1203  keyStart = "BostedChristy-p18@Pdg=";
1204  for(RgIMap::const_iterator it = entries.begin(); it != entries.end(); ++it)
1205  {
1206  const std::string& key = it->first;
1207  int pdg = 0;
1208  int A = 0;
1209  if (0 == key.compare(0, keyStart.size(), keyStart.c_str()))
1210  {
1211  pdg = atoi(key.c_str() + keyStart.size());
1212  A = pdg::IonPdgCodeToA(pdg);
1213  }
1214  if (0 != pdg && A != 0)
1215  {
1216  std::ostringstream key_ss ;
1217  key_ss << keyStart << pdg;
1218  RgKey rgkey = key_ss.str();
1219  double p18;
1220  GetParam( rgkey, p18) ;
1221  fMEC2009p18.insert(map<int,double>::value_type(A,p18));
1222  }
1223  }
1224 
1225 }
std::array< double, 7 > fMassRes
resonance mass
const int kPdgP33m1232_DeltaPP
Definition: PDGCodes.h:107
std::array< std::array< double, 4 >, 7 > fRescoefTp
tunable parameters from Ref.1, Table III for resonance
std::array< std::array< double, 5 >, 2 > fNRcoefTp
tunable parameters from Ref.1, Table III for nonres bkg
std::array< double, 9 > fEMCalpha
tunable parameters for EMC fit
const int kPdgF15m1680_NP
Definition: PDGCodes.h:135
#define pFATAL
Definition: Messenger.h:56
std::array< double, 8 > fMEC2009coef
tunable parameters for MEC2009 function
int GetParamVect(const std::string &comm_name, std::vector< T > &v, bool is_top_call=true) const
Handle to load vectors of parameters.
int IonPdgCodeToA(int pdgc)
Definition: PDGUtils.cxx:63
std::array< double, 3 > fEMCc
tunable parameters for EMC fit
double Mass(Resonance_t res)
resonance mass (GeV)
const int kPdgS11m1650_N0
Definition: PDGCodes.h:112
std::array< std::array< double, 3 >, 7 > fBRD
branching ratios of resonances for deterium fit
double Width(Resonance_t res)
resonance width (GeV)
const int kPdgF37m1950_DeltaM
Definition: PDGCodes.h:148
const int kPdgF37m1950_DeltaP
Definition: PDGCodes.h:150
const int kPdgS11m1535_N0
Definition: PDGCodes.h:108
virtual const Registry & GetConfig(void) const
Definition: Algorithm.cxx:246
std::array< double, 6 > fMECcoef
tunable parameters for Eqs.(20), (21) Ref.2
const int kPdgP33m1232_DeltaP
Definition: PDGCodes.h:106
const int kPdgP33m1232_DeltaM
Definition: PDGCodes.h:104
#define pALERT
Definition: Messenger.h:57
std::array< double, 13 > fAfitcoef
tunable parameters for nuclei fit
const int kPdgF37m1950_Delta0
Definition: PDGCodes.h:149
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:96
const RgIMap & GetItemMap(void) const
Definition: Registry.h:161
static constexpr double A
Definition: Units.h:74
void BranchingRatios(int, double &, double &) const
const int kPdgEta
Definition: PDGCodes.h:161
const int kPdgPi0
Definition: PDGCodes.h:160
int OrbitalAngularMom(Resonance_t res)
orbital angular momentum
const int kPdgF37m1950_DeltaPP
Definition: PDGCodes.h:151
std::array< std::array< double, 3 >, 7 > fBRp
branching ratios of resonances for proton fit
const int kPdgD13m1520_NP
Definition: PDGCodes.h:111
const int kPdgF15m1680_N0
Definition: PDGCodes.h:134
std::array< double, 7 > fWidthRes
resonance width
const int kPdgS11m1650_NP
Definition: PDGCodes.h:113
Resonance_t FromPdgCode(int pdgc)
PDG code -&gt; resonance id.
const int kPdgP33m1232_Delta0
Definition: PDGCodes.h:105
const int kPdgP11m1440_N0
Definition: PDGCodes.h:126
std::array< std::array< double, 5 >, 2 > fNRcoefTD
tunable parameters from Ref.1, Table IV for nonres bkg
const int kPdgD13m1520_N0
Definition: PDGCodes.h:110
const int kPdgP11m1440_NP
Definition: PDGCodes.h:127
static PDGLibrary * Instance(void)
Definition: PDGLibrary.cxx:68
bool fUseMEC
account for MEC contribution?
Singleton class to load &amp; serve a TDatabasePDG.
Definition: PDGLibrary.h:35
string RgKey
std::array< std::array< double, 4 >, 7 > fRescoefTD
tunable parameters from Ref.2, Table III for resonance
std::array< double, 6 > fNRcoefL
tunable parameters from Ref.1, Table III for nonres bkg
const int kPdgProton
Definition: PDGCodes.h:81
std::array< int, 7 > fAngRes
resonance angular momentum
TParticlePDG * Find(int pdgc, bool must_exist=true)
Definition: PDGLibrary.cxx:86
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 int kPdgS11m1535_NP
Definition: PDGCodes.h:109
const int kPdgTgtDeuterium
Definition: PDGCodes.h:201
std::array< std::array< double, 3 >, 7 > fRescoefL
tunable parameters from Ref.1, Table III for resonance
map< RgKey, RegistryItemI * > RgIMap
Definition: Registry.h:45
double BostedChristyEMPXSec::MEC2009 ( int  A,
double  Q2,
double  W 
) const
private

Definition at line 407 of file BostedChristyEMPXSec.cxx.

References fAM, fMEC2009coef, fMEC2009p18, genie::utils::kinematics::Q2(), and genie::utils::kinematics::W().

Referenced by XSec().

408 {
409  double F1 = 0.0;
410  double W2 = W*W;
411  double Mp = fAM;
412  double Mp2 = Mp*Mp;
413  if(W2<=0.0)
414  return F1;
415  double nu = (W2 - Mp2 + Q2)/2./Mp;
416  double x = Q2/2.0/Mp/nu;
417 
418  if(A<=2)
419  return F1;
420 
421  double p18;
422  for (const auto& kv : fMEC2009p18)
423  {
424  p18 = kv.second;
425  if (A<=kv.first)
426  break;
427  }
428 
429  F1 = fMEC2009coef[0]*TMath::Exp(-(W - fMEC2009coef[1])*(W - fMEC2009coef[1])/fMEC2009coef[2])/
430  TMath::Power(1. + TMath::Max(0.3,Q2)/fMEC2009coef[3],fMEC2009coef[4])*TMath::Power(nu, fMEC2009coef[5])*(1.0 +
431  p18*TMath::Power(A, fMEC2009coef[6] + fMEC2009coef[7]*x));
432 
433  if(F1<=1.0E-9)
434  F1 = 0.0;
435 
436  return F1;
437 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
std::array< double, 8 > fMEC2009coef
tunable parameters for MEC2009 function
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
static constexpr double A
Definition: Units.h:74
double BostedChristyEMPXSec::sigmaNR ( int  sf,
double  Q2,
double  W,
bool  isDeuterium = false 
) const
private

Definition at line 155 of file BostedChristyEMPXSec.cxx.

References fMP, fMpi0, fNRcoefL, fNRcoefTD, fNRcoefTp, genie::units::ub, and genie::utils::kinematics::W().

Referenced by FermiSmearingA(), FermiSmearingD(), and XSec().

156 {
157  const std::array<std::array<double, 5>, 2> &fNRcoefT = !isDeuterium?fNRcoefTp:fNRcoefTD;
158  if (isDeuterium)
159  sf=0;
160  double W2 = W*W;
161  double Mp = fMP;
162  double Mpi = fMpi0;
163  double Mp2 = Mp*Mp;
164 
165  double Wdif = W - (Mp + Mpi);
166 
167  double m0 = (sf==0) ? 0.125 : 4.2802; //Ref.1, Eqs.(22, 24)
168 
169  double Q20 = (sf==0) ? 0.05 : 0.125; //Ref.1, Eqs.(22, 24)
170 
171  double xpr = 1./(1.+(W2-(Mp+Mpi)*(Mp+Mpi))/(Q2+Q20)); // Ref.1, Eq.(22)
172 
173  double xsec = 0.0;
174 
175 
176  if (sf==0)
177  {
178 
179  for (int i=0;i<2;i++)
180  {
181  double h_nr = fNRcoefT[i][0]/TMath::Power(Q2+fNRcoefT[i][1], fNRcoefT[i][2]+fNRcoefT[i][3]*Q2+fNRcoefT[i][4]*Q2*Q2); // Ref.1, Eq. (21)
182  xsec += h_nr*TMath::Power(Wdif, 1.5+i);
183  }
184 
185  xsec *= xpr;
186  }
187  else
188  {
189  double xb = Q2/(Q2+W2-Mp2);
190  double t = TMath::Log(TMath::Log((Q2+m0)/0.330/0.330)/TMath::Log(m0/0.330/0.330)); // Ref.1, Eq. (24)
191  xsec += fNRcoefL[0]*TMath::Power(1.-xpr, fNRcoefL[2]+fNRcoefL[1]*t)/(1.-xb)/(Q2+Q20)
192  *TMath::Power(Q2/(Q2+Q20), fNRcoefL[3])*TMath::Power(xpr, fNRcoefL[4]+fNRcoefL[5]*t); // Ref.1, Eq. (23)
193  }
194 
195  return xsec*units::ub;
196 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
std::array< std::array< double, 5 >, 2 > fNRcoefTp
tunable parameters from Ref.1, Table III for nonres bkg
static constexpr double ub
Definition: Units.h:80
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
std::array< std::array< double, 5 >, 2 > fNRcoefTD
tunable parameters from Ref.1, Table IV for nonres bkg
std::array< double, 6 > fNRcoefL
tunable parameters from Ref.1, Table III for nonres bkg
double BostedChristyEMPXSec::sigmaR ( int  sf,
double  Q2,
double  W,
bool  isDeuterium = false 
) const
private

Definition at line 53 of file BostedChristyEMPXSec.cxx.

References genie::units::A, fAngRes, fBRD, fBRp, fMassRes, fMeta, fMP, fMpi0, fRescoefL, fRescoefTD, fRescoefTp, fWidthRes, genie::utils::kinematics::Q2(), genie::units::ub, and genie::utils::kinematics::W().

Referenced by FermiSmearingA(), FermiSmearingD(), and XSec().

54 {
55  const std::array<std::array<double, 3>, 7> &fBR = !isDeuterium?fBRp:fBRD;
56  const std::array<std::array<double, 4>, 7> &fRescoefT = !isDeuterium?fRescoefTp:fRescoefTD;
57  if (isDeuterium)
58  sf=0;
59 
60  double W2 = W*W;
61  // proton mass
62  double Mp = fMP;
63  // pion mass
64  double Mpi = fMpi0;
65  // eta-meson mass
66  double Meta = fMeta;
67 
68  double Mp2 = Mp*Mp;
69  double Mpi2 = Mpi*Mpi;
70  double Meta2 = Meta*Meta;
71 
72  // Calculate kinematics needed for threshold Relativistic B-W
73 
74  // Ref.1, Eq. (10)
75  double k = (W2 - Mp2)/2./Mp;
76  // Ref.1, Eq. (11)
77  double kcm = (W2 - Mp2)/2./W;
78  // mesons energy and momentim
79  double Epicm = (W2 + Mpi2 - Mp2)/2./W; // pion energy in CMS
80  double ppicm = TMath::Sqrt(TMath::Max(0.0, Epicm*Epicm - Mpi2)); // pion momentum in CMS
81  double Epi2cm = (W2 + 4*Mpi2 - Mp2)/2./W; // two pion energy in CMS
82  double ppi2cm = TMath::Sqrt(TMath::Max(0.0, Epi2cm*Epi2cm - 4*Mpi2)); // two pion energi n CMS
83  double Eetacm = (W2 + Meta2 - Mp2 )/2./W; // eta energy in CMS
84  double petacm = TMath::Sqrt(TMath::Max(0.0, Eetacm*Eetacm - Meta2)); // eta energy in CMS
85 
86  double xsec = 0.0;
87 
88  // going through seven resonances
89  for (int i=0;i<7;i++)
90  {
91  // resonance mass squared
92  double MassRes2 = fMassRes[i]*fMassRes[i];
93  // resonance damping parameter
94  double x0 = i!=0?0.215:!isDeuterium?0.14462:0.1446;
95  // Ref.1, Eq. (12)
96  double kr = (MassRes2-Mp2)/2./Mp;
97  // Ref.1, Eq. (13)
98  double kcmr = (MassRes2-Mp2)/2./fMassRes[i];
99 
100  // formulas analogous to the above with substitution W->MR_i
101  double Epicmr = (MassRes2 + Mpi2 - Mp2)/2./fMassRes[i];
102  double ppicmr = TMath::Sqrt(TMath::Max(0.0, Epicmr*Epicmr - Mpi2));
103  double Epi2cmr = (MassRes2 + 4.*Mpi2 - Mp2)/2./fMassRes[i];
104  double ppi2cmr = TMath::Sqrt(TMath::Max(0.0, Epi2cmr*Epi2cmr - 4.*Mpi2));
105  double Eetacmr = (MassRes2 + Meta2 - Mp2)/2./fMassRes[i];
106  double petacmr = TMath::Sqrt(TMath::Max(0.0, Eetacmr*Eetacmr - Meta2));
107 
108  // Calculate partial widths
109  // Ref.1, Eq. (15) for single pion
110  double pwid0 = fWidthRes[i]*TMath::Power(ppicm/ppicmr, 1.+2.*fAngRes[i])*
111  TMath::Power((ppicmr*ppicmr + x0*x0)/(ppicm*ppicm+x0*x0), fAngRes[i]); // 1-pion decay mode
112  // Ref.1, Eq. (16) for two pions
113  double pwid1 = 0;
114  if (!isDeuterium || (isDeuterium && i!=1))
115  pwid1 = W/fMassRes[i]*fWidthRes[i]*TMath::Power(ppi2cm/ppi2cmr, 4.+2.*fAngRes[i])*
116  TMath::Power((ppi2cmr*ppi2cmr + x0*x0)/(ppi2cm*ppi2cm + x0*x0), 2.+fAngRes[i]); // 2-pion decay mode
117  else
118  pwid1 = fWidthRes[i]*TMath::Power(petacm/petacmr, 1.+2.*fAngRes[i])*TMath::Power((ppi2cmr*ppi2cmr + x0*x0)/(ppi2cm*ppi2cm + x0*x0), fAngRes[i]);
119 
120 
121  double pwid2 = 0.; // eta decay mode
122  // Ref.1, Eq. (15) for eta
123  if(!isDeuterium && (i==1 || i==4))
124  pwid2 = fWidthRes[i]*TMath::Power(petacm/petacmr, 1.+2.*fAngRes[i])*TMath::Power((petacmr*petacmr + x0*x0)/(petacm*petacm + x0*x0), fAngRes[i]); // eta decay only for S11's
125 
126  // Ref.1, Eq. (17)
127  double pgam = fWidthRes[i]*(kcm/kcmr)*(kcm/kcmr)*(kcmr*kcmr+x0*x0)/(kcm*kcm+x0*x0);
128  // Ref.1, Eq. (14)
129  double width = fBR[i][0]*pwid0+fBR[i][1]*pwid1+fBR[i][2]*pwid2;
130 
131  // Begin resonance Q^2 dependence calculations
132  double A;
133 
134  if (sf==0)
135  // Ref.1, Eq. (18)
136  A = fRescoefT[i][0]*(1.+fRescoefT[i][1]*Q2/(1.+fRescoefT[i][2]*Q2))/TMath::Power(1.+Q2/0.91, fRescoefT[i][3]);
137  else
138  // Ref.1, Eq. (19)
139  A = fRescoefL[i][0]*Q2/(1.+fRescoefL[i][1]*Q2)*TMath::Exp(-1.*fRescoefL[i][2]*Q2);
140 
141 
142  // Ref.1, Eq. (9)
143  double BW = kr/k*kcmr/kcm/fWidthRes[i]*width*pgam/((W2 - MassRes2)*(W2 - MassRes2) + MassRes2*width*width);
144 
145  // Ref.1, Eq. (8) divided by W
146  xsec += BW*A*A;
147  }
148  // restore factor W in Ref.1, Eq. (8)
149  return W*xsec*units::ub;
150 }
std::array< double, 7 > fMassRes
resonance mass
std::array< std::array< double, 4 >, 7 > fRescoefTp
tunable parameters from Ref.1, Table III for resonance
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
std::array< std::array< double, 3 >, 7 > fBRD
branching ratios of resonances for deterium fit
static constexpr double ub
Definition: Units.h:80
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
static constexpr double A
Definition: Units.h:74
std::array< std::array< double, 3 >, 7 > fBRp
branching ratios of resonances for proton fit
std::array< double, 7 > fWidthRes
resonance width
std::array< std::array< double, 4 >, 7 > fRescoefTD
tunable parameters from Ref.2, Table III for resonance
std::array< int, 7 > fAngRes
resonance angular momentum
std::array< std::array< double, 3 >, 7 > fRescoefL
tunable parameters from Ref.1, Table III for resonance
bool BostedChristyEMPXSec::ValidKinematics ( const Interaction i) const
virtual

Is the input kinematical point a physically allowed one?

Reimplemented from genie::XSecAlgorithmI.

Definition at line 644 of file BostedChristyEMPXSec.cxx.

References fQ2max, fWmax, genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), genie::Kinematics::W(), and genie::utils::kinematics::W().

Referenced by XSec().

645 {
646  const Kinematics & kinematics = interaction -> Kine();
647  double W = kinematics.W();
648  double Q2 = kinematics.Q2();
649  if (W<fWmin || W>fWmax)
650  return false;
651  if (Q2<fQ2min || Q2>fQ2max)
652  return false;
653 
654  return true;
655 }
double W(bool selected=false) const
Definition: Kinematics.cxx:157
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
bool BostedChristyEMPXSec::ValidProcess ( const Interaction i) const
virtual

Can this cross section algorithm handle the input process?

Implements genie::XSecAlgorithmI.

Definition at line 619 of file BostedChristyEMPXSec.cxx.

References genie::Target::HitNucPdg(), genie::Interaction::InitState(), genie::pdg::IsChargedLepton(), genie::ProcessInfo::IsEM(), genie::pdg::IsNeutron(), genie::pdg::IsProton(), genie::ProcessInfo::IsResonant(), genie::kISkipProcessChk, genie::InitialState::ProbePdg(), genie::Interaction::ProcInfo(), and genie::InitialState::Tgt().

Referenced by XSec().

620 {
621  if(interaction->TestBit(kISkipProcessChk)) return true;
622 
623  const InitialState & init_state = interaction->InitState();
624  const ProcessInfo & proc_info = interaction->ProcInfo();
625 
626  if(!proc_info.IsResonant() ) return false;
627 
628 
629  int hitnuc = init_state.Tgt().HitNucPdg();
630  bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc));
631 
632  if (!is_pn) return false;
633 
634  int probe = init_state.ProbePdg();
635  bool is_em = proc_info.IsEM();
636 
637  bool l_em = (pdg::IsChargedLepton(probe) && is_em );
638 
639  if (!l_em) return false;
640 
641  return true;
642 }
bool IsResonant(void) const
Definition: ProcessInfo.cxx:99
int HitNucPdg(void) const
Definition: Target.cxx:304
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
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 BostedChristyEMPXSec::XSec ( const Interaction i,
KinePhaseSpace_t  k 
) const
virtual

Compute the cross section for the input interaction.

Implements genie::XSecAlgorithmI.

Definition at line 439 of file BostedChristyEMPXSec.cxx.

References genie::Target::A(), genie::units::A, fAfitcoef, FermiSmearingA(), FermiSmearingD(), FitEMC(), fKFTable, fMP, fNucRmvE, fPM, genie::utils::mec::J(), genie::utils::kinematics::Jacobian(), genie::constants::kAem, genie::constants::kAem2, genie::constants::kPi, genie::constants::kPi2, genie::kPSWQ2fE, genie::kRfHitNucRest, MEC2009(), genie::InitialState::ProbeE(), genie::Kinematics::Q2(), genie::utils::kinematics::Q2(), sigmaNR(), sigmaR(), genie::InitialState::Tgt(), ValidKinematics(), ValidProcess(), genie::Kinematics::W(), genie::utils::kinematics::W(), and genie::Target::Z().

441 {
442  if(! this -> ValidProcess (interaction) ) return 0.;
443  if(! this -> ValidKinematics (interaction) ) return 0.;
444 
445  // Get kinematical parameters
446  const Kinematics & kinematics = interaction -> Kine();
447  const InitialState & init_state = interaction -> InitState();
448  const Target & target = init_state.Tgt();
449  int A = target.A();
450  int Z = target.Z();
451  double E = init_state.ProbeE(kRfHitNucRest);
452  double W = kinematics.W();
453  double Q2 = kinematics.Q2();
454  double Wsq = W*W;
455  // Cross section for proton or neutron
456 
457  double Mp = fMP;
458  double Mp2 = Mp*Mp;
459  double MN = fPM;
460  double MN2 = MN*MN;
461 
462  double nu = (Wsq - MN2 + Q2)/2./MN;
463  double x = Q2/2./MN/nu;
464 
465  double sigmaT, sigmaL, F1p, R, W1;
466  // Cross section for proton or neutron
467  if (A<2 && Wsq>1.155)
468  {
469  double xb = Q2/(Wsq+Q2-Mp2);
470  sigmaT = sigmaR(0, Q2, W) + sigmaNR(0, Q2, W);
471  sigmaL = sigmaR(1, Q2, W) + sigmaNR(1, Q2, W);
472  F1p = sigmaT*(Wsq-Mp2)/8./kPi2/kAem;
473  R = sigmaL/sigmaT;
474  // If neutron, subtract proton from deuteron. Factor of two to
475  // convert from per nucleon to per deuteron
476  if(Z==0)
477  {
478  sigmaT = sigmaR(0, Q2, W, true) + sigmaNR(0, Q2, W, true);
479  double F1d = sigmaT*(Wsq-Mp2)/8./kPi2/kAem;
480  F1p = 2.*F1d - F1p;
481  }
482  W1 = F1p/MN;
483  }
484 
485  // For deuteron
486  if(A==2)
487  {
488  double Rd, F1c, F1d;
489  //get Fermi-smeared R from Erics proton fit
490  FermiSmearingD(Q2, W, F1c, R, sigmaT, sigmaL);
491  //get fit to F1 in deuteron, per nucleon
492  FermiSmearingD(Q2, W, F1d, Rd, sigmaT, sigmaL, true);
493  //convert to W1 per deuteron
494  W1 = F1d/MN*2.0;
495  }
496 
497  //For nuclei
498  if (A>2)
499  {
500  // Modifed to use Superscaling from Ref. 3
501  double Es, pF, kF;
502  for (const auto& kv : fNucRmvE)
503  {
504  Es = kv.second;
505  if (A<=kv.first)
506  break;
507  }
508  for (const auto& kv : fKFTable)
509  {
510  kF = kv.second;
511  if (A<=kv.first)
512  break;
513  }
514  // adjust pf to give right width based on kf
515  pF = 0.5*kF;
516  double F1d;
517  FermiSmearingA(Q2, W, pF, Es, F1p, F1d, sigmaT, sigmaL);
518  R = 0.;
519  if(sigmaT>0.)
520  R = sigmaL/sigmaT;
521  W1 = (2.*Z*F1d + (A - 2.*Z)*(2.*F1d - F1p))/MN;
522 
523  W1 *= (fAfitcoef[0] + x*(fAfitcoef[1] + x*(fAfitcoef[2] + x*(fAfitcoef[3] + x*(fAfitcoef[4] + x*fAfitcoef[5])))));
524 
525  if(W>0.)
526  W1 *= TMath::Power(fAfitcoef[6] + (fAfitcoef[7]*W + fAfitcoef[8]*Wsq)/(fAfitcoef[9] + fAfitcoef[10]*Q2),2);
527 
528  double F1M = MEC2009(A, Q2, W);
529 
530  W1 += F1M;
531  if(Wsq>0.)
532  R *= (fAfitcoef[11] + fAfitcoef[12]*A);
533  }
534 
535  double emcfac = FitEMC(x, A);
536 
537  W1 *= emcfac;
538 
539  double nu2 = nu*nu;
540  double Eprime = E - nu;
541  double sin2theta_2 = Q2/4/E/Eprime;
542  double cos2theta_2 = 1 - sin2theta_2;
543  double W2 = W1*(1 + R)/ (1+nu2/Q2);
544  double xsec = 4*Eprime*Eprime*kAem2/Q2/Q2*(2*W1*sin2theta_2 + W2*cos2theta_2); // d2xsec/dOmegadEprime
545  double jacobian = W*kPi/E/Eprime/MN;
546  xsec*= jacobian; // d2xsec/dOmegadEprime-> d2xsec/dWdQ2
547 
548  // The algorithm computes d^2xsec/dWdQ2
549  // Check whether variable tranformation is needed
550  if(kps!=kPSWQ2fE) {
551  double J = utils::kinematics::Jacobian(interaction,kPSWQ2fE,kps);
552  xsec *= J;
553  }
554 
555  return xsec;
556 
557 }
double W(bool selected=false) const
Definition: Kinematics.cxx:157
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int A(void) const
Definition: Target.h:70
double sigmaR(int, double, double, bool) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
void FermiSmearingA(double, double, double, double, double &, double &, double &, double &) const
double W(const Interaction *const i)
Definition: KineUtils.cxx:1101
std::array< double, 13 > fAfitcoef
tunable parameters for nuclei fit
static constexpr double A
Definition: Units.h:74
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
int Z(void) const
Definition: Target.h:68
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
double FitEMC(double, int) const
double MEC2009(int, double, double) const
double Q2(bool selected=false) const
Definition: Kinematics.cxx:125
const Target & Tgt(void) const
Definition: InitialState.h:66
void FermiSmearingD(double, double, double &, double &, double &, double &, bool) const
bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
double ProbeE(RefFrame_t rf) const
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
double sigmaNR(int, double, double, bool) const
Initial State information.
Definition: InitialState.h:48

Member Data Documentation

std::array<double, 13> genie::BostedChristyEMPXSec::fAfitcoef
private

tunable parameters for nuclei fit

Definition at line 100 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::BostedChristyEMPXSec::fAM
private

mass parameter

Definition at line 73 of file BostedChristyEMPXSec.h.

Referenced by FermiSmearingD(), LoadConfig(), and MEC2009().

std::array<int, 7> genie::BostedChristyEMPXSec::fAngRes
private

resonance angular momentum

Definition at line 85 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

std::array<std::array<double, 3>, 7> genie::BostedChristyEMPXSec::fBRD
private

branching ratios of resonances for deterium fit

Definition at line 83 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

std::array<std::array<double, 3>, 7> genie::BostedChristyEMPXSec::fBRp
private

branching ratios of resonances for proton fit

Definition at line 82 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

std::array<double, 9> genie::BostedChristyEMPXSec::fEMCalpha
private

tunable parameters for EMC fit

Definition at line 102 of file BostedChristyEMPXSec.h.

Referenced by FitEMC(), and LoadConfig().

std::array<double, 3> genie::BostedChristyEMPXSec::fEMCc
private

tunable parameters for EMC fit

Definition at line 103 of file BostedChristyEMPXSec.h.

Referenced by FitEMC(), and LoadConfig().

map<int, double> genie::BostedChristyEMPXSec::fKFTable
private

Definition at line 106 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and XSec().

std::array<double, 7> genie::BostedChristyEMPXSec::fMassRes
private

resonance mass

Definition at line 87 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

double genie::BostedChristyEMPXSec::fMD
private

deuterium mass

Definition at line 74 of file BostedChristyEMPXSec.h.

Referenced by FermiSmearingD(), and LoadConfig().

std::array<double, 8> genie::BostedChristyEMPXSec::fMEC2009coef
private

tunable parameters for MEC2009 function

Definition at line 99 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and MEC2009().

map<int, double> genie::BostedChristyEMPXSec::fMEC2009p18
private

Definition at line 105 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and MEC2009().

std::array<double, 6> genie::BostedChristyEMPXSec::fMECcoef
private

tunable parameters for Eqs.(20), (21) Ref.2

Definition at line 98 of file BostedChristyEMPXSec.h.

Referenced by FermiSmearingD(), and LoadConfig().

double genie::BostedChristyEMPXSec::fMeta
private

eta mass

Definition at line 76 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

double genie::BostedChristyEMPXSec::fMP
private

mass parameter

Definition at line 72 of file BostedChristyEMPXSec.h.

Referenced by FermiSmearingA(), FermiSmearingD(), LoadConfig(), sigmaNR(), sigmaR(), and XSec().

double genie::BostedChristyEMPXSec::fMpi0
private

pion mass

Definition at line 75 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), sigmaNR(), and sigmaR().

std::array<double, 6> genie::BostedChristyEMPXSec::fNRcoefL
private

tunable parameters from Ref.1, Table III for nonres bkg

Definition at line 97 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaNR().

std::array<std::array<double, 5>, 2> genie::BostedChristyEMPXSec::fNRcoefTD
private

tunable parameters from Ref.1, Table IV for nonres bkg

Definition at line 96 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaNR().

std::array<std::array<double, 5>, 2> genie::BostedChristyEMPXSec::fNRcoefTp
private

tunable parameters from Ref.1, Table III for nonres bkg

Definition at line 95 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaNR().

map<int, double> genie::BostedChristyEMPXSec::fNucRmvE
private

Definition at line 107 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and XSec().

double genie::BostedChristyEMPXSec::fPM
private

mass parameter

Definition at line 71 of file BostedChristyEMPXSec.h.

Referenced by FermiSmearingA(), LoadConfig(), and XSec().

double genie::BostedChristyEMPXSec::fQ2max
private

maximal Q2

Definition at line 80 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and ValidKinematics().

double genie::BostedChristyEMPXSec::fQ2min
private

minimal Q2

Definition at line 79 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig().

std::array<std::array<double, 3>, 7> genie::BostedChristyEMPXSec::fRescoefL
private

tunable parameters from Ref.1, Table III for resonance

Definition at line 93 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

std::array<std::array<double, 4>, 7> genie::BostedChristyEMPXSec::fRescoefTD
private

tunable parameters from Ref.2, Table III for resonance

Definition at line 92 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

std::array<std::array<double, 4>, 7> genie::BostedChristyEMPXSec::fRescoefTp
private

tunable parameters from Ref.1, Table III for resonance

Definition at line 91 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

bool genie::BostedChristyEMPXSec::fUseMEC
private

account for MEC contribution?

Definition at line 70 of file BostedChristyEMPXSec.h.

Referenced by FermiSmearingD(), and LoadConfig().

std::array<double, 7> genie::BostedChristyEMPXSec::fWidthRes
private

resonance width

Definition at line 89 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and sigmaR().

double genie::BostedChristyEMPXSec::fWmax
private

maximal W

Definition at line 78 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig(), and ValidKinematics().

double genie::BostedChristyEMPXSec::fWmin
private

minimal W

Definition at line 77 of file BostedChristyEMPXSec.h.

Referenced by LoadConfig().

const XSecIntegratorI* genie::BostedChristyEMPXSec::fXSecIntegrator
private

Definition at line 109 of file BostedChristyEMPXSec.h.

Referenced by Integral().


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