23 #include "Framework/Conventions/GBuild.h" 34 using namespace genie;
39 ROOT::Math::IBaseFunctionOneDim(),
64 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 65 LOG(
"GSLXSecFunc",
pDEBUG) <<
"xsec(Q2 = " << Q2 <<
") = " << xsec;
69 ROOT::Math::IBaseFunctionOneDim *
78 ROOT::Math::IBaseFunctionOneDim(),
102 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 103 LOG(
"GXSecFunc",
pDEBUG) <<
"xsec(y = " << y <<
") = " << xsec;
107 ROOT::Math::IBaseFunctionOneDim *
116 double DNuMass,
double scale) :
117 ROOT::Math::IBaseFunctionOneDim(),
141 double DNuEnergy = xin;
146 double E_nu = P4_nu->E();
149 double ETimesM = E_nu * M_target;
150 double EPlusM = E_nu + M_target;
152 double p_DNu = TMath::Sqrt(DNuEnergy*DNuEnergy - fDNuMass2);
153 double cos_theta_DNu = (DNuEnergy*(EPlusM) - ETimesM - 0.5*fDNuMass2) / (E_nu * p_DNu);
154 double theta_DNu = TMath::ACos(cos_theta_DNu);
155 TVector3 DNu_3vector = TVector3(0,0,0);
156 DNu_3vector.SetMagThetaPhi(p_DNu, theta_DNu, 0.);
157 TLorentzVector P4_DNu = TLorentzVector(DNu_3vector, DNuEnergy);
160 TVector3 target_3vector = P4_nu->Vect() - DNu_3vector;
161 double E_target = E_nu + M_target - DNuEnergy;
162 TLorentzVector P4_target = TLorentzVector(target_3vector , E_target);
164 kinematics->
SetQ2(2.*M_target*(E_target-M_target));
168 #ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ 169 LOG(
"GSLXSecFunc",
pDEBUG) <<
"xsec(DNuEnergy = " << DNuEnergy <<
") = " << xsec;
174 ROOT::Math::IBaseFunctionOneDim *
185 const double M2 = M * M;
188 const double A = M2 + 2.*M*
E;
189 const double B = (M+
E) * (E*M + 0.5*fDNuMass2);
190 const double C = E*E *(M2 + fDNuMass2) + E*M*fDNuMass2 + 0.25*fDNuMass2*fDNuMass2;
191 const double D = sqrt(B*B - A*C);
193 Range1D_t DNuEnergy((B - D)/A, (B + D)/A);
200 ROOT::Math::IBaseFunctionMultiDim(),
230 ROOT::Math::IBaseFunctionMultiDim *
239 ROOT::Math::IBaseFunctionMultiDim(),
269 ROOT::Math::IBaseFunctionMultiDim *
278 ROOT::Math::IBaseFunctionMultiDim(),
307 ROOT::Math::IBaseFunctionMultiDim *
316 ROOT::Math::IBaseFunctionMultiDim(),
350 ROOT::Math::IBaseFunctionMultiDim *
359 ROOT::Math::IBaseFunctionMultiDim(),
392 ROOT::Math::IBaseFunctionMultiDim *
401 ROOT::Math::IBaseFunctionMultiDim(),
440 ROOT::Math::IBaseFunctionMultiDim *
449 ROOT::Math::IBaseFunctionOneDim(),
477 ROOT::Math::IBaseFunctionOneDim *
486 ROOT::Math::IBaseFunctionOneDim(),
514 ROOT::Math::IBaseFunctionOneDim *
523 ROOT::Math::IBaseFunctionOneDim(),
551 ROOT::Math::IBaseFunctionOneDim *
560 ROOT::Math::IBaseFunctionOneDim(),
588 ROOT::Math::IBaseFunctionOneDim *
600 ROOT::Math::IBaseFunctionMultiDim(),
626 double E_nu = P4_nu->E();
629 double theta_l = xin[1];
630 double phi_l = xin[2];
631 double theta_pi = xin[3];
632 double phi_pi = xin[4];
634 double E_pi= E_nu-E_l;
636 double y = E_pi/E_nu;
652 double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
653 TVector3 lepton_3vector = TVector3(0,0,0);
654 lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
655 TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
657 double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
658 TVector3 pion_3vector = TVector3(0,0,0);
659 pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
660 TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
662 double Q2 = -(*P4_nu-P4_lep).Mag2();
668 if ( x < xlim.min || x > xlim.
max ) {
680 if (xsec>0 &&
flip) {
688 ROOT::Math::IBaseFunctionMultiDim *
701 ROOT::Math::IBaseFunctionMultiDim(),
725 double E_nu = P4_nu->E();
728 double theta_l = xin[1];
729 double phi_l = xin[2];
730 double theta_pi = xin[3];
731 double phi_pi = xin[4];
733 double E_pi= E_nu-E_l;
735 double y = E_pi/E_nu;
750 double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
751 TVector3 lepton_3vector = TVector3(0,0,0);
752 lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
753 TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
755 double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
756 TVector3 pion_3vector = TVector3(0,0,0);
757 pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
758 TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
760 double Q2 = -(*P4_nu-P4_lep).Mag2();
766 if ( x < xlim.min || x > xlim.
max ) {
783 ROOT::Math::IBaseFunctionMultiDim *
798 ROOT::Math::IBaseFunctionMultiDim(),
824 double E_nu = P4_nu->E();
827 double theta_l = xin[1];
829 double theta_pi = xin[2];
830 double phi_pi = xin[3];
832 double sin_theta_l = TMath::Sin(theta_l);
833 double sin_theta_pi = TMath::Sin(theta_pi);
835 double E_pi= E_nu-E_l;
837 double y = E_pi/E_nu;
852 double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
853 TVector3 lepton_3vector = TVector3(0,0,0);
854 lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
855 TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
857 double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
858 TVector3 pion_3vector = TVector3(0,0,0);
859 pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
860 TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
862 double Q2 = -(*P4_nu-P4_lep).Mag2();
868 if ( x < xlim.min || x > xlim.
max ) {
884 ROOT::Math::IBaseFunctionMultiDim *
901 ROOT::Math::IBaseFunctionMultiDim(),
927 double E_nu = P4_nu->E();
931 double theta_l = xin[0];
932 double phi_l = xin[1];
933 double theta_pi = xin[2];
936 double sin_theta_l = TMath::Sin(theta_l);
937 double sin_theta_pi = TMath::Sin(theta_pi);
939 double E_pi= E_nu-E_l;
941 double y = E_pi/E_nu;
956 double p_l = TMath::Sqrt(E_l*E_l - m_l*m_l);
957 TVector3 lepton_3vector = TVector3(0,0,0);
958 lepton_3vector.SetMagThetaPhi(p_l,theta_l,phi_l);
959 TLorentzVector P4_lep = TLorentzVector(lepton_3vector , E_l );
961 double p_pi = TMath::Sqrt(E_pi*E_pi - m_pi*m_pi);
962 TVector3 pion_3vector = TVector3(0,0,0);
963 pion_3vector.SetMagThetaPhi(p_pi,theta_pi,phi_pi);
964 TLorentzVector P4_pion = TLorentzVector(pion_3vector , E_pi);
966 double Q2 = -(*P4_nu-P4_lep).Mag2();
972 if ( x < xlim.min || x > xlim.
max ) {
1003 string gsl_nd_integrator_type,
double gsl_relative_tolerance,
1004 unsigned int max_n_calls) :
1005 ROOT::Math::IBaseFunctionOneDim(),
1009 fGSLIntegratorType(gsl_nd_integrator_type),
1010 fGSLRelTol(gsl_relative_tolerance),
1011 fGSLMaxCalls(max_n_calls)
1015 integrator.SetRelTolerance(gsl_relative_tolerance);
1032 LOG(
"GSLXSecFunc",
pINFO) <<
"dXSec_dElep_AR >> "<<
func->
NDim()<<
"d integral done. (Elep = " <<Elep<<
" , dxsec/dElep = "<<xsec <<
")";
1043 const ROOT::Math::IBaseFunctionMultiDim *
fn,
1044 bool * ifLog,
double *
mins,
double *
maxes) :
1064 double * toEval =
new double[this->
NDim()];
1066 for (
unsigned int i = 0 ; i < this->
NDim() ; i++ )
1078 double val = (*fFn)(toEval);
const XSecAlgorithmI * fModel
unsigned int NDim(void) const
Cross Section Calculation Interface.
const XSecAlgorithmI * fModel
double DoEval(const double *xin) const
unsigned int NDim(void) const
const KPhaseSpace & PhaseSpace(void) const
static const double kNapierConst
const Interaction * fInteraction
unsigned int NDim(void) const
const XSecAlgorithmI * fModel
double DoEval(const double *xin) const
const Interaction * fInteraction
double DoEval(double xin) const
bool IsWeakCC(void) const
const Interaction * fInteraction
Range1D_t IntegrationRange(void) const
THE MAIN GENIE PROJECT NAMESPACE
d2XSec_dQ2dydt_E(const XSecAlgorithmI *m, const Interaction *i)
double DoEval(double xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
double DoEval(const double *xin) const
static const double kNucleonMass
unsigned int fGSLMaxCalls
double Q2(const Interaction *const i)
ROOT::Math::IntegrationMultiDim::Type IntegrationNDimTypeFromString(string type)
unsigned int NDim(void) const
void SetQ2(double Q2, bool selected=false)
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
d3Xsec_dOmegaldThetapi * Clone(void) const
Kinematics * KinePtr(void) const
dXSec_dQ2_E(const XSecAlgorithmI *m, const Interaction *i, double scale=1.)
virtual double XSec(const Interaction *i, KinePhaseSpace_t k=kPSfE) const =0
Compute the cross section for the input interaction.
double DoEval(double xin) const
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
A simple [min,max] interval for doubles.
static const double kPi0Mass
unsigned int NDim(void) const
const Interaction * fInteraction
double DoEval(double xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
Generated/set kinematical variables for an event.
double DoEval(const double *xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const Interaction * fInteraction
const XSecAlgorithmI * fModel
const XSecAlgorithmI * fModel
d5XSecAR(const XSecAlgorithmI *m, const Interaction *i)
unsigned int NDim(void) const
double DoEval(const double *xin) const
unsigned int NDim(void) const
const Interaction * fInteraction
const Interaction * fInteraction
double DoEval(double xin) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const Interaction * fInteraction
void UpdateWYFromXQ2(const Interaction *in)
const XSecAlgorithmI * fModel
double DoEval(double xin) const
double DoEval(const double *xin) const
d3Xsec_dOmegaldThetapi(const XSecAlgorithmI *m, const Interaction *i)
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
double DoEval(double xin) const
d2XSec_dxdy_Ex(const XSecAlgorithmI *m, const Interaction *i, double x)
unsigned int NDim(void) const
const XSecAlgorithmI * fModel
dXSec_dEDNu_E(const XSecAlgorithmI *m, const Interaction *i, double DNuMass, double scale=1.)
Summary information for an interaction.
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
unsigned int NDim(void) const
double DoEval(double xin) const
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
void SetFSLeptonP4(const TLorentzVector &p4)
static constexpr double cm2
~d4Xsec_dEldThetaldOmegapi()
const XSecAlgorithmI * fModel
d3XSec_dxdydt_E(const XSecAlgorithmI *m, const Interaction *i)
unsigned int NDim(void) const
dXSec_Log_Wrapper(const ROOT::Math::IBaseFunctionMultiDim *fn, bool *ifLog, double *min, double *maxes)
const XSecAlgorithmI * fModel
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
d2XSec_dWdQ2_EQ2(const XSecAlgorithmI *m, const Interaction *i, double Q2)
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
void Sett(double t, bool selected=false)
double DoEval(const double *xin) const
const XSecAlgorithmI * fModel
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
const XSecAlgorithmI * fModel
const Interaction * fInteraction
~d5Xsec_dEldOmegaldOmegapi()
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
dXSec_dy_E(const XSecAlgorithmI *m, const Interaction *i)
static const double kASmallNum
void SetE_lep(double E_lepton) const
unsigned int NDim(void) const
d2XSec_dQ2dy_E(const XSecAlgorithmI *m, const Interaction *i)
const Interaction * fInteraction
double DoEval(const double *xin) const
unsigned int NDim(void) const
TParticlePDG * FSPrimLepton(void) const
final state primary lepton
const Interaction * fInteraction
const XSecAlgorithmI * fModel
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
void SetFactor(double factor)
unsigned int NDim(void) const
const Interaction * fInteraction
const XSecAlgorithmI * fModel
bool IsDeepInelastic(void) const
void Setx(double x, bool selected=false)
d2XSec_dxdy_Ey(const XSecAlgorithmI *m, const Interaction *i, double x)
const XSecAlgorithmI * fModel
void UpdateWQ2FromXY(const Interaction *in)
void SetW(double W, bool selected=false)
TLorentzVector * HitNucP4Ptr(void) const
const Interaction * fInteraction
const Interaction * fInteraction
d2XSec_dxdy_E(const XSecAlgorithmI *m, const Interaction *i)
static const double kPionMass
string fGSLIntegratorType
const genie::utils::gsl::d3Xsec_dOmegaldThetapi * func
double DoEval(const double *xin) const
const XSecAlgorithmI * fModel
const XSecAlgorithmI * fModel
d5Xsec_dEldOmegaldOmegapi(const XSecAlgorithmI *m, const Interaction *i)
void Sety(double y, bool selected=false)
void UpdateXFromQ2Y(const Interaction *in)
unsigned int NDim(void) const
d4Xsec_dEldThetaldOmegapi(const XSecAlgorithmI *m, const Interaction *i)
const Interaction * fInteraction
const Interaction * fInteraction
~d2XSec_dlog10xdlog10Q2_E()
unsigned int NDim(void) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
d2XSec_dlog10xdlog10Q2_E(const XSecAlgorithmI *m, const Interaction *i, double scale=1.)
d2XSec_dWdQ2_E(const XSecAlgorithmI *m, const Interaction *i)
void SetHadSystP4(const TLorentzVector &p4)
const Interaction * fInteraction
double DoEval(const double *xin) const
InitialState * InitStatePtr(void) const
Range1D_t XLim(void) const
x limits
~d3Xsec_dOmegaldThetapi()
const InitialState & InitState(void) const
const ProcessInfo & ProcInfo(void) const
unsigned int NDim(void) const
double DoEval(const double *xin) const
const ROOT::Math::IBaseFunctionMultiDim * fFn
const Target & Tgt(void) const
ROOT::Math::IBaseFunctionMultiDim * Clone(void) const
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
unsigned int NDim(void) const
d2XSec_dWdQ2_EW(const XSecAlgorithmI *m, const Interaction *i, double W)
const Interaction * fInteraction
ROOT::Math::IBaseFunctionOneDim * Clone(void) const
ROOT::Math::IntegratorMultiDim integrator
bool IsDarkMatterDeepInelastic(void) const
dXSec_dElep_AR * Clone(void) const
double ProbeE(RefFrame_t rf) const
const XSecAlgorithmI * fModel
TLorentzVector * GetProbeP4(RefFrame_t rf=kRfHitNucRest) const
const XSecAlgorithmI * fModel
unsigned int NDim(void) const