ReinDFRPXSec.cxx
Go to the documentation of this file.
1 //____________________________________________________________________________
2 /*
3  Copyright (c) 2003-2019, The GENIE Collaboration
6
7  Author: Costas Andreopoulos <costas.andreopoulos \at stfc.ac.uk>
8  University of Liverpool & STFC Rutherford Appleton Lab - February 15, 2008
9
10  For the class documentation see the corresponding header file.
11
12  Important revisions after version 2.0.0 :
13  @ Feb 15, 2009 - CA
14  This class was first added in version 2.5.1.
15 */
16 //____________________________________________________________________________
17
18 #include <TMath.h>
19
21
32
33 using namespace genie;
34 using namespace genie::constants;
35 using namespace genie::controls;
36 using namespace genie::utils;
37
38 //____________________________________________________________________________
40 XSecAlgorithmI("genie::ReinDFRPXSec")
41 {
42
43 }
44 //____________________________________________________________________________
46 XSecAlgorithmI("genie::ReinDFRPXSec", config)
47 {
48
49 }
50 //____________________________________________________________________________
52 {
53
54 }
55 //____________________________________________________________________________
57  const Interaction * interaction, KinePhaseSpace_t kps) const
58 {
59  if(! this -> ValidProcess (interaction) ) return 0.;
60  if(! this -> ValidKinematics (interaction) ) return 0.;
61
62  const Kinematics & kinematics = interaction -> Kine();
63  const InitialState & init_state = interaction -> InitState();
64  const Target & target = init_state.Tgt();
65
66  bool isCC = interaction->ProcInfo().IsWeakCC();
67  double E = init_state.ProbeE(kRfHitNucRest); // neutrino energy
68  double x = kinematics.x(); // bjorken x
69  double y = kinematics.y(); // inelasticity y
70  double t = kinematics.t(); // (magnitude of) square of four-momentum xferred to proton
71  double M = target.HitNucMass(); //
72  double Q2 = 2.*x*y*M*E; // momentum transfer Q2>0
73  double Gf = kGF2 * M/(16*kPi3); // GF/pi/etc factor
74  double fp = 0.93 * kPionMass; // pion decay constant (cc)
75  double fp2 = TMath::Power(fp,2.);
76  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
77  double b = fBeta;
78  double ma2 = TMath::Power(fMa,2);
79  double propg = TMath::Power(ma2/(ma2+Q2),2.); // propagator term
80  double sTot = utils::hadxs::TotalPionNucleonXSec(Epi, isCC); // pi+N total cross section; CC process always produces a charged pion
81  double sTot2 = TMath::Power(sTot,2.);
82  double tFac = TMath::Exp(-b*t);
83
84 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
85  LOG("ReinDFR", pDEBUG)
86  << "E = " << E << ", x = " << x << ", y = " << y << ", Q2 = " << Q2;
87  LOG("ReinDFR", pDEBUG)
88  << "Epi = " << Epi << ", s^{piN}_{tot} = " << sTot;
89  LOG("ReinDFR", pDEBUG)
90  << "b = " << b << ", t = [" << tmin << ", " << tmax << "]";
91 #endif
92
93  // fixme: WARNING: don't leave this in here!!!
94  // only needed for comparison to Rein's paper
95 // double W2 = M*M + 2 * M * y * E - Q2;
96 // if (W2 < 4)
97 // return 0;
98
99  //----- compute d^2sigma/dxdydt
100  double xsec = Gf*E*fp2*(1-y)*propg*sTot2*tFac;
101
102  // NC XS is half of CC
103  if (!isCC)
104  xsec *= 0.5;
105
106  //----- Check whether variable tranformation is needed
107  if(kps!=kPSxytfE) {
108  double J = utils::kinematics::Jacobian(interaction,kPSxytfE,kps);
109 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__
110  LOG("ReinDFR", pDEBUG)
111  << "Jacobian for transformation to: "
112  << KinePhaseSpace::AsString(kps) << ", J = " << J;
113 #endif
114  xsec *= J;
115  }
116
117  //----- if requested return the free nucleon xsec even for input nuclear tgt
118  if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec;
119
120  //----- number of scattering centers in the target
121  int nucpdgc = target.HitNucPdg();
122  int NNucl = (pdg::IsProton(nucpdgc)) ? target.Z() : target.N();
123  xsec *= NNucl;
124
125  return xsec;
126 }
127 //____________________________________________________________________________
128 double ReinDFRPXSec::Integral(const Interaction * interaction) const
129 {
130  double xsec = fXSecIntegrator->Integrate(this,interaction);
131  return xsec;
132 }
133 //____________________________________________________________________________
134 bool ReinDFRPXSec::ValidProcess(const Interaction * interaction) const
135 {
136  //expect only free protons as target
137  const InitialState & init_state = interaction -> InitState();
138  const Target & target = init_state.Tgt();
139  if(target.A() > 1 || target.Z() != 1)
140  return false;
141
142  if(interaction->TestBit(kISkipProcessChk)) return true;
143
144  if(interaction->ProcInfo().IsDiffractive()) return true;
145  return false;
146 }
147 //____________________________________________________________________________
149 {
150  Algorithm::Configure(config);
152 }
153 //____________________________________________________________________________
155 {
156  Algorithm::Configure(config);
158 }
159 //____________________________________________________________________________
161 {
162  GetParam( "DFR-Ma", fMa ) ;
163  GetParam( "DFR-Beta",fBeta ) ;
164
166  dynamic_cast<const XSecIntegratorI *> (this->SubAlg("XSec-Integrator") );
168 }
169 //____________________________________________________________________________
170
Cross Section Calculation Interface.
Basic constants.
bool IsWeakCC(void) const
const XML_Char * target
Definition: expat.h:268
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
THE MAIN GENIE PROJECT NAMESPACE
Definition: GeneratorBase.h:8
Cross Section Integrator Interface.
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
int HitNucPdg(void) const
Definition: Target.cxx:321
int A(void) const
Definition: Target.h:71
double HitNucMass(void) const
Definition: Target.cxx:250
double Integral(const Interaction *i) const
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double x(bool selected=false) const
Definition: Kinematics.cxx:109
Definition: config.py:1
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:122
Summary information for an interaction.
Definition: Interaction.h:56
static const double kPi3
Definition: Constants.h:40
virtual bool ValidKinematics(const Interaction *i) const
Is the input kinematical point a physically allowed one?
bool IsProton(int pdgc)
Definition: PDGUtils.cxx:299
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Float_t E
Definition: plot.C:20
static string AsString(KinePhaseSpace_t kps)
A Neutrino Interaction Target. Is a transparent encapsulation of quite different physical systems suc...
Definition: Target.h:41
virtual void Configure(const Registry &config)
Definition: Algorithm.cxx:70
const XSecIntegratorI * fXSecIntegrator
Definition: ReinDFRPXSec.h:56
int Z(void) const
Definition: Target.h:69
Misc GENIE control constants.
Double_t xsec[nknots]
Definition: testXsec.C:47
int N(void) const
Definition: Target.h:70
static const double kPionMass
Definition: Constants.h:74
A registry. Provides the container for algorithm configuration parameters.
Definition: Registry.h:66
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.
virtual double Integrate(const XSecAlgorithmI *model, const Interaction *interaction) const =0
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
Definition: KineUtils.cxx:128
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:180
FILE * fp
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:67
static const double kGF2
Definition: Constants.h:60
void Configure(const Registry &config)
void kinematics()
Definition: kinematics.C:10
double ProbeE(RefFrame_t rf) const
Root of GENIE utility namespaces.
const UInt_t kISkipProcessChk
if set, skip process validity checks
Definition: Interaction.h:47
double TotalPionNucleonXSec(double Epion, bool isChargedPion=true)