GENIEGenerator
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
ReinDFRPXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2024, The GENIE Collaboration
4  For the full text of the license visit http://copyright.genie-mc.org
5 
6  Costas Andreopoulos <c.andreopoulos \at cern.ch>
7  University of Liverpool
8 */
9 //____________________________________________________________________________
10 
11 #include <TMath.h>
12 
15 #include "Framework/Conventions/GBuild.h"
25 
26 using namespace genie;
27 using namespace genie::constants;
28 using namespace genie::controls;
29 using namespace genie::utils;
30 
31 //____________________________________________________________________________
33 XSecAlgorithmI("genie::ReinDFRPXSec")
34 {
35 
36 }
37 //____________________________________________________________________________
38 ReinDFRPXSec::ReinDFRPXSec(const std::string & config) :
39 XSecAlgorithmI("genie::ReinDFRPXSec", config)
40 {
41 
42 }
43 //____________________________________________________________________________
45 {
46 
47 }
48 //____________________________________________________________________________
50  const Interaction * interaction, KinePhaseSpace_t kps) const
51 {
52  if(! this -> ValidProcess (interaction) ) return 0.;
53  if(! this -> ValidKinematics (interaction) ) return 0.;
54 
55  const Kinematics & kinematics = interaction -> Kine();
56  const InitialState & init_state = interaction -> InitState();
57  const Target & target = init_state.Tgt();
58 
59  bool isCC = interaction->ProcInfo().IsWeakCC();
60  double E = init_state.ProbeE(kRfHitNucRest); // neutrino energy
61  double x = kinematics.x(); // bjorken x
62  double y = kinematics.y(); // inelasticity y
63  double t = kinematics.t(); // (magnitude of) square of four-momentum xferred to proton
64  double M = target.HitNucMass(); //
65  double Q2 = 2.*x*y*M*E; // momentum transfer Q2>0
66  double Gf = kGF2 * M/(16*kPi3); // GF/pi/etc factor
67  double fp = 0.93 * kPionMass; // pion decay constant (cc)
68  double fp2 = TMath::Power(fp,2.);
69  double Epi = y*E - t/(2*M); // pion energy. note we use - instead of + like Rein's paper b/c our t is magnitude only
70  double b = fBeta;
71  double ma2 = TMath::Power(fMa,2);
72  double propg = TMath::Power(ma2/(ma2+Q2),2.); // propagator term
73  double sTot = utils::hadxs::TotalPionNucleonXSec(Epi, isCC); // pi+N total cross section; CC process always produces a charged pion
74  double sTot2 = TMath::Power(sTot,2.);
75  double tFac = TMath::Exp(-b*t);
76 
77 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
78  LOG("ReinDFR", pDEBUG)
79  << "E = " << E << ", x = " << x << ", y = " << y << ", Q2 = " << Q2;
80  LOG("ReinDFR", pDEBUG)
81  << "Epi = " << Epi << ", s^{piN}_{tot} = " << sTot;
82  LOG("ReinDFR", pDEBUG)
83  << "b = " << b << ", t = [" << tmin << ", " << tmax << "]";
84 #endif
85 
86  // fixme: WARNING: don't leave this in here!!!
87  // only needed for comparison to Rein's paper
88 // double W2 = M*M + 2 * M * y * E - Q2;
89 // if (W2 < 4)
90 // return 0;
91 
92  //----- compute d^2sigma/dxdydt
93  double xsec = Gf*E*fp2*(1-y)*propg*sTot2*tFac;
94 
95  // NC XS is half of CC
96  if (!isCC)
97  xsec *= 0.5;
98 
99  //----- Check whether variable tranformation is needed
100  if(kps!=kPSxytfE) {
101  double J = utils::kinematics::Jacobian(interaction,kPSxytfE,kps);
102 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
103  LOG("ReinDFR", pDEBUG)
104  << "Jacobian for transformation to: "
105  << KinePhaseSpace::AsString(kps) << ", J = " << J;
106 #endif
107  xsec *= J;
108  }
109 
110  //----- if requested return the free nucleon xsec even for input nuclear tgt
111  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
112 
113  //----- number of scattering centers in the target
114  int nucpdgc = target.HitNucPdg();
115  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
116  xsec *= NNucl;
117 
118  return xsec;
119 }
120 //____________________________________________________________________________
121 double ReinDFRPXSec::Integral(const Interaction * interaction) const
122 {
123  double xsec = fXSecIntegrator->Integrate(this,interaction);
124  return xsec;
125 }
126 //____________________________________________________________________________
127 bool ReinDFRPXSec::ValidProcess(const Interaction * interaction) const
128 {
129  //expect only free protons as target
130  const InitialState & init_state = interaction -> InitState();
131  const Target & target = init_state.Tgt();
132  if(target.A() > 1 || target.Z() != 1)
133  return false;
134 
135  if(interaction->TestBit(kISkipProcessChk)) return true;
136 
137  if(interaction->ProcInfo().IsDiffractive()) return true;
138  return false;
139 }
140 //____________________________________________________________________________
141 void ReinDFRPXSec::Configure(const Registry & config)
142 {
143  Algorithm::Configure(config);
144  this->LoadConfig();
145 }
146 //____________________________________________________________________________
147 void ReinDFRPXSec::Configure(string config)
148 {
149  Algorithm::Configure(config);
150  this->LoadConfig();
151 }
152 //____________________________________________________________________________
154 {
155  GetParam( "DFR-Ma", fMa ) ;
156  GetParam( "DFR-Beta",fBeta ) ;
157 
159  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator") );
160  assert(fXSecIntegrator);
161 }
162 //____________________________________________________________________________
Cross Section Calculation Interface.
bool IsWeakCC(void) const
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:147
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:1077
int HitNucPdg(void) const
Definition: Target.cxx:304
int A(void) const
Definition: Target.h:70
double HitNucMass(void) const
Definition: Target.cxx:233
double Integral(const Interaction *i) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:39
double x(bool selected=false) const
Definition: Kinematics.cxx:99
bool IsDiffractive(void) const
enum genie::EKinePhaseSpace KinePhaseSpace_t
double fBeta
b in dsig{piN}/dt = dsig0{piN}/dt * exp(-b(t-tmin)), b ~ 0.333 (nucleon_size)^2
Definition: ReinDFRPXSec.h:54
double y(bool selected=false) const
Definition: Kinematics.cxx:112
static constexpr double b
Definition: Units.h:78
Summary information for an interaction.
Definition: Interaction.h:56
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
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
static string AsString(KinePhaseSpace_t kps)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:40
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:62
const XSecIntegratorI * fXSecIntegrator
Definition: ReinDFRPXSec.h:56
int Z(void) const
Definition: Target.h:68
int N(void) const
Definition: Target.h:69
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:65
const UInt_t kIAssumeFreeNucleon
Definition: Interaction.h:49
double XSec(const Interaction *i, KinePhaseSpace_t k) const
Compute the cross section for the input interaction.
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:130
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double fMa
axial mass
Definition: ReinDFRPXSec.h:53
double t(bool selected=false) const
Definition: Kinematics.cxx:170
bool ValidProcess(const Interaction *i) const
Can this cross section algorithm handle the input process?
bool GetParam(const RgKey &name, T &p, bool is_top_call=true) const
const Target & Tgt(void) const
Definition: InitialState.h:66
void Configure(const Registry &config)
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
double ProbeE(RefFrame_t rf) const
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)
Definition: HadXSUtils.cxx:51
Initial State information.
Definition: InitialState.h:48
#define pDEBUG
Definition: Messenger.h:63
const Algorithm * SubAlg(const RgKey &registry_key) const
Definition: Algorithm.cxx:345