genie::utils::kinematics Namespace Reference

Kinematical utilities. More...

Namespaces
namespace	electromagnetic

Functions
double	PhaseSpaceVolume (const Interaction *const i, KinePhaseSpace_t ps)
double	Jacobian (const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
bool	TransformMatched (KinePhaseSpace_t ia, KinePhaseSpace_t ib, KinePhaseSpace_t a, KinePhaseSpace_t b, bool &fwd)
Range1D_t	InelWLim (double Ev, double M, double ml)
Range1D_t	InelQ2Lim_W (double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Range1D_t	Inelq2Lim_W (double Ev, double M, double ml, double W, double q2min_cut=-1 *controls::kMinQ2Limit)
Range1D_t	InelQ2Lim (double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Range1D_t	Inelq2Lim (double Ev, double M, double ml, double q2min_cut=-1 *controls::kMinQ2Limit)
Range1D_t	InelXLim (double Ev, double M, double ml)
Range1D_t	InelYLim (double Ev, double M, double ml)
Range1D_t	InelYLim_X (double Ev, double M, double ml, double x)
Range1D_t	CohW2Lim (double Mn, double m_produced, double mlep, double Ev, double Q2)
Range1D_t	CohNuLim (double W2min, double W2max, double Q2, double Mn, double xsi)
Range1D_t	CohYLim (double Mn, double m_produced, double mlep, double Ev, double Q2, double xsi)
Range1D_t	CohYLim (double EvL, double ml)
Range1D_t	CohXLim (void)
Range1D_t	CohQ2Lim (double Mn, double m_produced, double mlep, double Ev)
Range1D_t	Cohq2Lim (double Mn, double m_produced, double mlep, double Ev)
Range1D_t	CEvNSQ2Lim (double Ev)
Range1D_t	DarkWLim (double Ev, double M, double ml)
Range1D_t	DarkQ2Lim_W (double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Range1D_t	Darkq2Lim_W (double Ev, double M, double ml, double W, double q2min_cut=-1 *controls::kMinQ2Limit)
Range1D_t	DarkQ2Lim (double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Range1D_t	Darkq2Lim (double Ev, double M, double ml, double q2min_cut=-1 *controls::kMinQ2Limit)
Range1D_t	DarkXLim (double Ev, double M, double ml)
Range1D_t	DarkYLim (double Ev, double M, double ml)
Range1D_t	DarkYLim_X (double Ev, double M, double ml, double x)
double	CohW2Min (double Mn, double m_produced)
double	QD2toQ2 (double QD2)
double	Q2toQD2 (double Q2)
void	WQ2toXY (double Ev, double M, double W, double Q2, double &x, double &y)
void	XYtoWQ2 (double Ev, double M, double &W, double &Q2, double x, double y)
void	XQ2toWY (double Ev, double M, double &W, double Q2, double x, double &y)
double	XYtoW (double Ev, double M, double x, double y)
double	XYtoQ2 (double Ev, double M, double x, double y)
double	Q2YtoX (double Ev, double M, double Q2, double y)
void	UpdateWQ2FromXY (const Interaction *in)
void	UpdateXYFromWQ2 (const Interaction *in)
void	UpdateWYFromXQ2 (const Interaction *in)
void	UpdateXFromQ2Y (const Interaction *in)
void	ApplyCutsToKineLimits (Range1D_t &r, double min, double max)
double	Q2 (const Interaction *const i)
double	W (const Interaction *const i)
bool	IsAboveCharmThreshold (double x, double Q2, double M, double mc)
double	SlowRescalingVar (double x, double Q2, double M, double mc)
double	RESImportanceSamplingEnvelope (double *x, double *par)
double	DISImportanceSamplingEnvelope (double *x, double *par)
double	COHImportanceSamplingEnvelope (double *x, double *par)

Detailed Description

Kinematical utilities.

Author

Costas Andreopoulos <c.andreopoulos \at cern.ch> University of Liverpool

Changes required to implement the GENIE Boosted Dark Matter module were installed by Josh Berger (Univ. of Wisconsin)

Created:\n November 26, 2004

License:\n Copyright (c) 2003-2025, The GENIE Collaboration

For the full text of the license visit http://copyright.genie-mc.org

Function Documentation

ApplyCutsToKineLimits()

void genie::utils::kinematics::ApplyCutsToKineLimits	(	Range1D_t &	r,
		double	min,
		double	max )

Definition at line 1271 of file KineUtils.cxx.

{
  // if the min,max are within the existing limits, the cut can be applied
  // by narrowing down the xisting limits
  if ( utils::math::IsWithinLimits(min_cut, range ) ) range.min = min_cut;
  if ( utils::math::IsWithinLimits(max_cut, range ) ) range.max = max_cut;
 
  // if the min-cut is above the existing max-limit or
  // if the max-cut is below the existing min-limit then
  // the range should be invalidated
 
  if (min_cut > range.max || max_cut < range.min) {
 
    range.min = 0;
    range.max = 0;
  }
}

References genie::utils::math::IsWithinLimits(), genie::Range1D_t::max, and genie::Range1D_t::min.

CEvNSQ2Lim()

Range1D_t genie::utils::kinematics::CEvNSQ2Lim

(

double

Ev

)

Definition at line 886 of file KineUtils.cxx.

{
  Range1D_t Q2(controls::kAVerySmallNum, 4*Ev*Ev);
  return Q2;
}

References genie::controls::kAVerySmallNum, and Q2().

Referenced by genie::KPhaseSpace::Q2Lim().

COHImportanceSamplingEnvelope()

double genie::utils::kinematics::COHImportanceSamplingEnvelope	(	double *	x,
		double *	par )

Definition at line 1466 of file KineUtils.cxx.

{
  //-- inputs
  double xb = x[0];   // x
  double yb = x[1];   // y
 
  //-- parameters
  double xsmax = 3*par[0]; // safety factor * max cross section in (x,y)
  double Ev    = par[1]; // neutrino energy;
 
  if(yb<0.|| yb>1.) return 0.;
  if(xb<0.|| xb>1.) return 0.;
 
  if(Ev<1) return xsmax;
  if(xb/Ev<1E-4 && yb>0.95) return 5*xsmax;
 
  double func = 0;
  double xp   = 0.1;
  double yp   = (Ev>2.5) ? 2.5/Ev : 1;
 
  if(xb>xp) {
    double xs0=0;
    if(yb<yp) {
      xs0 = xsmax;
    } else {
      xs0 = xsmax - (yb-yp)*xsmax;
    }
    double d = TMath::Power( (xb-xp)/0.075, 2);
    func = xs0/(1 + d);
  } else {
    if(yb>yp) {
       func = xsmax - (yb-yp)*xsmax;
    } else {
       func = xsmax;
    }
  }
  return func;
}

References func().

Referenced by genie::COHKinematicsGenerator::LoadConfig().

CohNuLim()

Range1D_t genie::utils::kinematics::CohNuLim	(	double	W2min,
		double	W2max,
		double	Q2,
		double	Mn,
		double	xsi )

Definition at line 823 of file KineUtils.cxx.

{
  Range1D_t nul;
  nul.min = -1;
  nul.max = -1;
 
  double nu_min = (W2min + Q2 - Mn * Mn) / (2.0 * Mn);
  double nu_max = (W2max + Q2 - Mn * Mn) / (2.0 * Mn);
  double xsiQ = xsi * TMath::Sqrt(Q2);
 
  nul.min = (xsiQ > nu_min) ? xsiQ : nu_min;
  nul.max = nu_max;
 
  return nul;
}

References genie::Range1D_t::max, genie::Range1D_t::min, and Q2().

Referenced by CohYLim().

CohQ2Lim()

Range1D_t genie::utils::kinematics::CohQ2Lim	(	double	Mn,
		double	m_produced,
		double	mlep,
		double	Ev )

Definition at line 743 of file KineUtils.cxx.

{
  // The expressions for Q^2 min appears in PRD 74, 054007 (2006) by
  // Kartavtsev, Paschos, and Gounaris
 
  // That expression is specified for the case of the pion.
  // GENIE has also Coherent interaction with Single gamma production and Rho
  // so that formula will be adapted for a generic m_produced
 
  Range1D_t Q2;
  Q2.min = 0.0;
  Q2.max = std::numeric_limits<double>::max();  // Value must be overriden in user options
 
  double Mn2 = Mn * Mn;
  double mlep2 = mlep * mlep;
  double s = Mn2 + 2.0 * Mn * Ev;
  double W2min = CohW2Min(Mn, m_produced);
 
  // Looks like Q2min = A * B - C, where A, B, and C are complicated
  double a = 1.0;
  double b = mlep2 / s;
  double c = W2min / s;
  double lambda = a * a + b * b + c * c - 2.0 * a * b - 2.0 * a * c - 2.0 * b * c;
  if (lambda > 0) {
    double A = (s - Mn * Mn) / 2.0;
    double B = 1 - TMath::Sqrt(lambda);
    double C = 0.5 * (W2min + mlep2 - Mn2 * (W2min - mlep2) / s );
    if (A * B - C < 0) {
      SLOG("KineLimits", pERROR)
        << "Q2 kinematic limits calculation failed for CohQ2Lim. "
        << "Assuming Q2min = 0.0";
    }
    Q2.min = TMath::Max(0., A * B - C);
  } else {
    SLOG("KineLimits", pERROR)
      << "Q2 kinematic limits calculation failed for CohQ2Lim. "
      << "Assuming Q2min = 0.0";
  }
 
  return Q2;
}

References a, CohW2Min(), pERROR, Q2(), and SLOG.

Referenced by Cohq2Lim(), and genie::KPhaseSpace::Q2Lim().

Cohq2Lim()

Range1D_t genie::utils::kinematics::Cohq2Lim	(	double	Mn,
		double	m_produced,
		double	mlep,
		double	Ev )

Definition at line 785 of file KineUtils.cxx.

{
  Range1D_t Q2 = utils::kinematics::CohQ2Lim(Mn, m_produced, mlep, Ev);
  Range1D_t q2;
  q2.min = - Q2.max;
  q2.max = - Q2.min;
  return q2;
}

References CohQ2Lim(), genie::Range1D_t::max, genie::Range1D_t::min, and Q2().

CohW2Lim()

Range1D_t genie::utils::kinematics::CohW2Lim	(	double	Mn,
		double	m_produced,
		double	mlep,
		double	Ev,
		double	Q2 )

Definition at line 794 of file KineUtils.cxx.

{
  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
  // Kartavtsev, Paschos, and Gounaris
 
  // That expression is specified for the case of the pion.
  // GENIE has also Coherent interaction with Single gamma production and Rho
  // so that formula will be adapted for a generic m_produced
 
  Range1D_t W2l;
  W2l.min = -1;
  W2l.max = -1;
 
  double s = Mn * Mn + 2.0 * Mn * Ev;
  double Mnterm = 1 - Mn * Mn / s;
  double Mlterm = 1 - mlep * mlep / s;
  // Here T1, T2 are generically "term 1" and "term 2" in a long expression
  double T1 = 0.25 * s * s * Mnterm * Mnterm * Mlterm;
  double T2 = Q2 - (0.5 * s * Mnterm) + (0.5 * mlep * mlep * Mnterm);
 
  W2l.min = CohW2Min(Mn, m_produced);
  W2l.max = (T1 - T2 * T2 ) *
    (1.0 / Mnterm) *
    (1.0 / (Q2 + mlep * mlep));
 
  return W2l;
}

References CohW2Min(), genie::Range1D_t::max, genie::Range1D_t::min, and Q2().

Referenced by CohYLim().

CohW2Min()

double genie::utils::kinematics::CohW2Min	(	double	Mn,
		double	m_produced )

Definition at line 874 of file KineUtils.cxx.

{
  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
  // Kartavtsev, Paschos, and Gounaris
 
  // That expression is specified for the case of the pion.
  // GENIE has also Coherent interaction with Single gamma production and Rho
  // so that formula will be adapted for a generic m_produced
 
  return (Mn + m_produced) * (Mn + m_produced);
}

Referenced by CohQ2Lim(), and CohW2Lim().

CohXLim()

Range1D_t genie::utils::kinematics::CohXLim

(

void

)

Definition at line 735 of file KineUtils.cxx.

{
// Computes x limits for coherent v interactions
 
  Range1D_t x(controls::kASmallNum, 1.-controls::kASmallNum);
  return x;
}

References genie::controls::kASmallNum.

Referenced by genie::KPhaseSpace::XLim().

CohYLim() [1/2]

Range1D_t genie::utils::kinematics::CohYLim	(	double	EvL,
		double	ml )

Definition at line 865 of file KineUtils.cxx.

{
// Computes y limits for coherent v interactions
 
  Range1D_t y(kPionMass/EvL + controls::kASmallNum,
              1.-ml/EvL - controls::kASmallNum);
  return y;
}

References genie::controls::kASmallNum, and genie::constants::kPionMass.

CohYLim() [2/2]

Range1D_t genie::utils::kinematics::CohYLim	(	double	Mn,
		double	m_produced,
		double	mlep,
		double	Ev,
		double	Q2,
		double	xsi )

Definition at line 840 of file KineUtils.cxx.

{
  Range1D_t ylim;
  ylim.min = -1;
  ylim.max = -1;
 
  Range1D_t W2lim = genie::utils::kinematics::CohW2Lim(Mn, m_produced, mlep, Ev, Q2);
  if (W2lim.min > W2lim.max) {
    LOG("KineLimits", pDEBUG)
      << "Kinematically forbidden region in CohYLim. W2min = " << W2lim.min
      << "; W2max =" << W2lim.max;
    LOG("KineLimits", pDEBUG)
      << "  Mn = " << Mn << "; m_had_sys = " << m_produced << "; mlep = "
      << mlep << "; Ev = " << Ev << "; Q2 = " << Q2;
    return ylim;
  }
  Range1D_t nulim = genie::utils::kinematics::CohNuLim(W2lim.min, W2lim.max,
      Q2, Mn, xsi);
  ylim.min = nulim.min / Ev;
  ylim.max = nulim.max / Ev;
 
  return ylim;
}

References CohNuLim(), CohW2Lim(), LOG, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, and Q2().

Referenced by genie::KPhaseSpace::YLim(), genie::KPhaseSpace::YLim(), and genie::KPhaseSpace::YLim_X().

DarkQ2Lim()

Range1D_t genie::utils::kinematics::DarkQ2Lim	(	double	Ev,
		double	M,
		double	ml,
		double	Q2min_cut = controls::kMinQ2Limit )

Definition at line 973 of file KineUtils.cxx.

{
// Computes Q2 (>0) limits irrespective of W for inelastic v interactions
 
  Range1D_t Q2;
  Q2.min = -1;
  Q2.max = -1;
 
  Range1D_t W  = utils::kinematics::DarkWLim(Ev,M,ml);
  if(W.min<0) return Q2;
 
  Q2 = utils::kinematics::DarkQ2Lim_W(Ev,M,ml,W.min,Q2min_cut);
  return Q2;
}

References DarkQ2Lim_W(), DarkWLim(), Q2(), and W().

Referenced by Darkq2Lim(), and genie::KPhaseSpace::Q2Lim().

Darkq2Lim()

Range1D_t genie::utils::kinematics::Darkq2Lim	(	double	Ev,
		double	M,
		double	ml,
		double	q2min_cut = -1*controls::kMinQ2Limit )

Definition at line 989 of file KineUtils.cxx.

{
// Computes Q2 (>0) limits irrespective of W for inelastic v interactions
 
  Range1D_t Q2 = utils::kinematics::DarkQ2Lim(Ev,M,ml,-1.*q2min_cut);
  Range1D_t q2;
  q2.min = - Q2.max;
  q2.max = - Q2.min;
  return q2;
}

References DarkQ2Lim(), genie::Range1D_t::max, genie::Range1D_t::min, and Q2().

DarkQ2Lim_W()

Range1D_t genie::utils::kinematics::DarkQ2Lim_W	(	double	Ev,
		double	M,
		double	ml,
		double	W,
		double	Q2min_cut = controls::kMinQ2Limit )

Definition at line 915 of file KineUtils.cxx.

{
// Computes Q2 limits (>0) @ the input W for inelastic v interactions
 
  Range1D_t Q2;
  Q2.min = -1;
  Q2.max = -1;
 
  double M2  = TMath::Power(M,  2.);
  double ml2 = TMath::Power(ml, 2.);
  double W2  = TMath::Power(W,  2.);
  double s   = M2 + 2*M*Ev + ml2;
  assert(s > 0);
  double sqs = TMath::Sqrt(s);
  double E1CM = (s + ml2 - M2) / (2.*sqs);
  double p1CM = TMath::Max(0., E1CM*E1CM - ml2);
  p1CM = TMath::Sqrt(p1CM);
  double E3CM = (s + ml2 - W2) / (2.*sqs);
  double p3CM = TMath::Max(0., E3CM*E3CM - ml2);
  p3CM = TMath::Sqrt(p3CM);
 
  SLOG("KineLimits", pDEBUG) << "s  = " << s;
  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
  SLOG("KineLimits", pDEBUG) << "M = " << M;
  SLOG("KineLimits", pDEBUG) << "W = " << W;
  SLOG("KineLimits", pDEBUG) << "E1_CM = " << E1CM;
  SLOG("KineLimits", pDEBUG) << "p1_CM = " << p1CM;
  SLOG("KineLimits", pDEBUG) << "E3_CM = " << E3CM;
  SLOG("KineLimits", pDEBUG) << "p3_CM = " << p3CM;
 
  Q2.min = TMath::Power(p3CM - p1CM,2) - TMath::Power((W2 - M2) / (2.*sqs),2);
  Q2.max = TMath::Power(p3CM + p1CM,2) - TMath::Power((W2 - M2) / (2.*sqs),2);
 
  SLOG("KineLimits", pDEBUG) << "Nominal Q^2 limits: " << Q2.min << " , " << Q2.max;
  // guard against overflows
  Q2.max = TMath::Max(0., Q2.max);
  Q2.min = TMath::Max(0., Q2.min);
 
  // limit the minimum Q2
  if(Q2.min < Q2min_cut) {Q2.min = Q2min_cut;     }
  if(Q2.max < Q2.min   ) {Q2.min = -1; Q2.max = -1;}
 
  return Q2;
}

References pDEBUG, Q2(), SLOG, and W().

Referenced by DarkQ2Lim(), Darkq2Lim_W(), DarkXLim(), DarkYLim(), genie::KPhaseSpace::Q2Lim(), and genie::KPhaseSpace::Q2Lim_W().

Darkq2Lim_W()

Range1D_t genie::utils::kinematics::Darkq2Lim_W	(	double	Ev,
		double	M,
		double	ml,
		double	W,
		double	q2min_cut = -1*controls::kMinQ2Limit )

Definition at line 961 of file KineUtils.cxx.

{
// Computes q2 (<0) limits @ the input W for inelastic v interactions
 
  Range1D_t Q2 = utils::kinematics::DarkQ2Lim_W(Ev,M,ml,W,-1.*q2min_cut);
  Range1D_t q2;
  q2.min = - Q2.max;
  q2.max = - Q2.min;
  return q2;
}

References DarkQ2Lim_W(), genie::Range1D_t::max, genie::Range1D_t::min, Q2(), and W().

DarkWLim()

Range1D_t genie::utils::kinematics::DarkWLim	(	double	Ev,
		double	M,
		double	ml )

Definition at line 892 of file KineUtils.cxx.

{
// Computes W limits for inelastic v interactions
//
  double M2 = TMath::Power(M,2);
  double ml2 = TMath::Power(ml,2);
  double s  = M2 + 2*M*Ev + ml2;
  assert (s>0);
 
  Range1D_t W;
  W.min  = kNeutronMass + kPhotontest;
//  W.min  = kNeutronMass + kPionMass;
  W.max  = TMath::Sqrt(s) - ml;
  if(W.max<=W.min) {
    W.min = -1;
    W.max = -1;
    return W;
  }
  W.min  += controls::kASmallNum;
  W.max  -= controls::kASmallNum;
  return W;
}

References genie::controls::kASmallNum, genie::constants::kNeutronMass, genie::constants::kPhotontest, and W().

Referenced by DarkQ2Lim(), DarkXLim(), DarkYLim(), DarkYLim_X(), and genie::KPhaseSpace::WLim().

DarkXLim()

Range1D_t genie::utils::kinematics::DarkXLim	(	double	Ev,
		double	M,
		double	ml )

Definition at line 1001 of file KineUtils.cxx.

{
// Computes Bjorken x limits for inelastic interactions
// For the dark matter case, it is relatively straightforward
  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
  double Wmin = Wl.min;
  double W2min = Wmin*Wmin;
  SLOG("KineLimits", pDEBUG) << "W^2_min = " << W2min;
  Range1D_t Q2l = utils::kinematics::DarkQ2Lim_W(Ev, M, ml, Wmin);
  SLOG("KineLimits", pDEBUG) << "Q^2 range : " << Q2l.min << " , " << Q2l.max;
  double M2 = M*M;
  Range1D_t x;
  x.min = Q2l.min / (Q2l.min + W2min - M2);
  x.max = Q2l.max / (Q2l.max + W2min - M2);
 
  SLOG("KineLimits", pDEBUG) << "x  = [" << x.min << ", " << x.max << "]";
  return x;
}

References DarkQ2Lim_W(), DarkWLim(), genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, and SLOG.

Referenced by genie::KPhaseSpace::XLim().

DarkYLim()

Range1D_t genie::utils::kinematics::DarkYLim	(	double	Ev,
		double	M,
		double	ml )

Definition at line 1021 of file KineUtils.cxx.

{
  // For dark inelastic scattering, can compute exactly and is much simpler
  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
  double Wmin = Wl.min;
  double W2min = Wmin*Wmin;
  Range1D_t Q2l = utils::kinematics::DarkQ2Lim_W(Ev, M, ml, Wmin);
  double M2 = M*M;
  Range1D_t y;
  y.min = (Q2l.min + W2min - M2) / (2*Ev*M);
  y.max = (Q2l.max + W2min - M2) / (2*Ev*M);
 
  SLOG("KineLimits", pDEBUG) << "y  = [" << y.min << ", " << y.max << "]";
  return y;
}

References DarkQ2Lim_W(), DarkWLim(), genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, and SLOG.

Referenced by genie::KPhaseSpace::YLim().

DarkYLim_X()

Range1D_t genie::utils::kinematics::DarkYLim_X	(	double	Ev,
		double	M,
		double	ml,
		double	x )

Definition at line 1037 of file KineUtils.cxx.

{
  // Computes y limits @ the input x for inelastic interactions
  // We hit y_min when W = W_min and y_max when Q^2 = Q_min^2 or Q^2 = Q_max^2
 
  Range1D_t y;
  y.min = -1;
  y.max = -1;
 
  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
  double Wmin = Wl.min;
  double W2min = Wmin*Wmin;
  double M2 = M*M;
  double ml2 = ml*ml;
  y.min = (W2min - M2) / (1.-x) / (2*Ev*M);
  y.max = 2.* M * x *(Ev*Ev - ml2) / Ev / (2. * M * Ev * x + M2 * x * x + ml2);
 
  return y;
}

References DarkWLim(), genie::Range1D_t::max, and genie::Range1D_t::min.

Referenced by genie::KPhaseSpace::YLim_X().

DISImportanceSamplingEnvelope()

double genie::utils::kinematics::DISImportanceSamplingEnvelope	(	double *	x,
		double *	par )

Definition at line 1434 of file KineUtils.cxx.

{
  //-- inputs
  double xb = x[0];   // scaling variable x brorken
  //double y  = x[1];   // inelasticity y
 
  //-- parameters
  double xpeak = par[0]; // x peak
  double xmin  = par[1]; // min x
  double xmax  = 1.;     // max x
  double xsmax = par[2]; // safety factor * max cross section in (x,y)
 
  double func = 0;
 
  if(xb < xpeak/20.) {
       //low-x falling edge
       double plateau_edge = xpeak/20.;
       double slope = xsmax/(xmin - plateau_edge);
       func = xsmax - slope * (xb - plateau_edge);
  } else if (xb > 2*xpeak) {
       //high-x falling edge
       double plateau_edge = 2*xpeak;
       double slope = xsmax/(xmax - plateau_edge);
       func = xsmax - slope * (xb - plateau_edge);
  } else {
       // plateau
       func = xsmax;
  }
  return func;
}

References func().

InelQ2Lim()

Range1D_t genie::utils::kinematics::InelQ2Lim	(	double	Ev,
		double	M,
		double	ml,
		double	Q2min_cut = controls::kMinQ2Limit )

Definition at line 429 of file KineUtils.cxx.

{
// Computes Q2 (>0) limits irrespective of W for inelastic v interactions
 
  Range1D_t Q2;
  Q2.min = -1;
  Q2.max = -1;
 
  Range1D_t W  = utils::kinematics::InelWLim(Ev,M,ml);
  if(W.min<0) return Q2;
 
  Q2 = utils::kinematics::InelQ2Lim_W(Ev,M,ml,W.min,Q2min_cut);
  return Q2;
}

References InelQ2Lim_W(), InelWLim(), Q2(), and W().

Referenced by Inelq2Lim(), and genie::KPhaseSpace::Q2Lim().

Inelq2Lim()

Range1D_t genie::utils::kinematics::Inelq2Lim	(	double	Ev,
		double	M,
		double	ml,
		double	q2min_cut = -1*controls::kMinQ2Limit )

Definition at line 445 of file KineUtils.cxx.

{
// Computes Q2 (>0) limits irrespective of W for inelastic v interactions
 
  Range1D_t Q2 = utils::kinematics::InelQ2Lim(Ev,M,ml,-1.*q2min_cut);
  Range1D_t q2;
  q2.min = - Q2.max;
  q2.max = - Q2.min;
  return q2;
}

References InelQ2Lim(), genie::Range1D_t::max, genie::Range1D_t::min, and Q2().

InelQ2Lim_W()

Range1D_t genie::utils::kinematics::InelQ2Lim_W	(	double	Ev,
		double	M,
		double	ml,
		double	W,
		double	Q2min_cut = controls::kMinQ2Limit )

Definition at line 379 of file KineUtils.cxx.

{
// Computes Q2 limits (>0) @ the input W for inelastic v interactions
 
  Range1D_t Q2;
  Q2.min = -1;
  Q2.max = -1;
 
  double M2  = TMath::Power(M,  2.);
  double ml2 = TMath::Power(ml, 2.);
  double W2  = TMath::Power(W,  2.);
  double s   = M2 + 2*M*Ev;
 
  SLOG("KineLimits", pDEBUG) << "s  = " << s;
  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
  assert (s>0);
 
  double auxC = 0.5*(s-M2)/s;
  double aux1 = s + ml2 - W2;
  double aux2 = aux1*aux1 - 4*s*ml2;
 
  (aux2 < 0) ? ( aux2 = 0 ) : ( aux2 = TMath::Sqrt(aux2) );
 
  Q2.max = -ml2 + auxC * (aux1 + aux2); // => 0
  Q2.min = -ml2 + auxC * (aux1 - aux2); // => 0
 
  // guard against overflows
  Q2.max = TMath::Max(0., Q2.max);
  Q2.min = TMath::Max(0., Q2.min);
 
  // limit the minimum Q2
  if(Q2.min < Q2min_cut) {Q2.min = Q2min_cut;     }
  if(Q2.max < Q2.min   ) {Q2.min = -1; Q2.max = -1;}
 
  return Q2;
}

References pDEBUG, Q2(), SLOG, and W().

Referenced by InelQ2Lim(), Inelq2Lim_W(), genie::KPhaseSpace::Q2Lim(), genie::KPhaseSpace::Q2Lim_W(), and genie::EmpiricalMECPXSec2015::XSec().

Inelq2Lim_W()

Range1D_t genie::utils::kinematics::Inelq2Lim_W	(	double	Ev,
		double	M,
		double	ml,
		double	W,
		double	q2min_cut = -1*controls::kMinQ2Limit )

Definition at line 417 of file KineUtils.cxx.

{
// Computes q2 (<0) limits @ the input W for inelastic v interactions
 
  Range1D_t Q2 = utils::kinematics::InelQ2Lim_W(Ev,M,ml,W,-1.*q2min_cut);
  Range1D_t q2;
  q2.min = - Q2.max;
  q2.max = - Q2.min;
  return q2;
}

References InelQ2Lim_W(), genie::Range1D_t::max, genie::Range1D_t::min, Q2(), and W().

InelWLim()

Range1D_t genie::utils::kinematics::InelWLim	(	double	Ev,
		double	M,
		double	ml )

Definition at line 358 of file KineUtils.cxx.

{
// Computes W limits for inelastic v interactions
//
  double M2 = TMath::Power(M,2);
  double s  = M2 + 2*M*Ev;
  assert (s>0);
 
  Range1D_t W;
  W.min  = kNeutronMass + kPhotontest;
  W.max  = TMath::Sqrt(s) - ml;
  if(W.max<=W.min) {
    W.min = -1;
    W.max = -1;
    return W;
  }
  W.min  += controls::kASmallNum;
  W.max  -= controls::kASmallNum;
  return W;
}

References genie::controls::kASmallNum, genie::constants::kNeutronMass, genie::constants::kPhotontest, and W().

Referenced by InelQ2Lim(), genie::KPhaseSpace::WLim(), and genie::EmpiricalMECPXSec2015::XSec().

InelXLim()

Range1D_t genie::utils::kinematics::InelXLim	(	double	Ev,
		double	M,
		double	ml )

Definition at line 457 of file KineUtils.cxx.

{
// Computes Bjorken x limits for inelastic v interactions
 
  double M2  = TMath::Power(M, 2.);
  double ml2 = TMath::Power(ml,2.);
  double s   = M2 + 2*M*Ev;
 
  SLOG("KineLimits", pDEBUG) << "s  = " << s;
  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
  assert (s>M2);
 
  Range1D_t x;
  x.min = ml2/(s-M2) + controls::kAVerySmallNum;
  x.max = 1.         - controls::kAVerySmallNum;
 
  return x;
}

References genie::controls::kAVerySmallNum, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, and SLOG.

Referenced by InelYLim(), and genie::KPhaseSpace::XLim().

InelYLim()

Range1D_t genie::utils::kinematics::InelYLim	(	double	Ev,
		double	M,
		double	ml )

Definition at line 477 of file KineUtils.cxx.

{
// Computes y limits for inelastic v interactions
 
  Range1D_t y;
  y.min =  999;
  y.max = -999;
 
  Range1D_t xl = kinematics::InelXLim(Ev,M,ml);
  assert(xl.min>0 && xl.max>0);
 
  const unsigned int N=100;
  const double logxmin = TMath::Log10(xl.min);
  const double logxmax = TMath::Log10(xl.max);
  const double dlogx   = (logxmax-logxmin) / (double)(N-1);
 
  for(unsigned int i=0; i<N; i++) {
    double x = TMath::Power(10, logxmin + i*dlogx);
 
    Range1D_t y_x = kinematics::InelYLim_X(Ev,M,ml,x);
    if(y_x.min>=0 && y_x.min<=1) y.min = TMath::Min(y.min, y_x.min);
    if(y_x.max>=0 && y_x.max<=1) y.max = TMath::Max(y.max, y_x.max);
  }
 
  if(y.max >= 0 && y.max <= 1 && y.min >= 0 && y.min <= 1) {
    y.min = TMath::Max(y.min,     controls::kAVerySmallNum);
    y.max = TMath::Min(y.max, 1 - controls::kAVerySmallNum);
  } else {
    y.min = -1;
    y.max = -1;
  }
  SLOG("KineLimits", pDEBUG) << "y  = [" << y.min << ", " << y.max << "]";
  return y;
}

References InelXLim(), InelYLim_X(), genie::controls::kAVerySmallNum, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, and SLOG.

Referenced by genie::KPhaseSpace::YLim().

InelYLim_X()

Range1D_t genie::utils::kinematics::InelYLim_X	(	double	Ev,
		double	M,
		double	ml,
		double	x )

Definition at line 512 of file KineUtils.cxx.

{
// Computes y limits @ the input x for inelastic v interactions
 
  Range1D_t y;
  y.min = -1;
  y.max = -1;
 
  double Ev2 = TMath::Power(Ev,2);
  double ml2 = TMath::Power(ml,2);
 
  SLOG("KineLimits", pDEBUG) << "x  = " << x;
  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
 
  assert (Ev>0);
  assert (x>0&&x<1);
 
  double a = 0.5 * ml2/(M*Ev*x);
  double b = ml2/Ev2;
  double c = 1 + 0.5*x*M/Ev;
  double d = TMath::Max(0., TMath::Power(1-a,2.) - b);
 
  double A = 0.5 * (1-a-0.5*b)/c;
  double B = 0.5 * TMath::Sqrt(d)/c;
 
  y.min = TMath::Max(0., A-B) + controls::kAVerySmallNum;
  y.max = TMath::Min(1., A+B) - controls::kAVerySmallNum;
 
  return y;
}

References a, genie::controls::kAVerySmallNum, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, and SLOG.

Referenced by InelYLim(), and genie::KPhaseSpace::YLim_X().

IsAboveCharmThreshold()

bool genie::utils::kinematics::IsAboveCharmThreshold	(	double	x,
		double	Q2,
		double	M,
		double	mc )

Definition at line 1238 of file KineUtils.cxx.

{
// x : scaling variable (plain or modified)
// Q2: momentum transfer
// M : hit nucleon "mass" (nucleon can be off the mass shell)
// mc: charm mass
//
  double M2   = TMath::Power(M,2);
  double v    = 0.5*Q2/(M*x);
  double W2   = TMath::Max(0., M2+2*M*v-Q2);
  double W    = TMath::Sqrt(W2);
  double Wmin = M + kLightestChmHad;
  double xc   = utils::kinematics::SlowRescalingVar(x,Q2,M,mc);
 
  if(xc>=1 || W<=Wmin) return false;
  else                 return true;
}

References genie::constants::kLightestChmHad, Q2(), SlowRescalingVar(), and W().

Referenced by genie::QPMDISStrucFuncBase::CalcPDFs(), and genie::QPMDMDISStrucFuncBase::CalcPDFs().

Jacobian()

double genie::utils::kinematics::Jacobian	(	const Interaction *const	i,
		KinePhaseSpace_t	f,
		KinePhaseSpace_t	t )

Definition at line 130 of file KineUtils.cxx.

{
// Returns the Jacobian for a kinematical transformation:
// from_ps (x,y,z,...) -> to_ps (u,v,w,...).
//
// Note on the convention used here:
// The differential cross-section d^n sigma / dxdydz..., is transformed to
// d^n sigma / dudvdw... using the Jacobian, computed in this function, as:
// as d^n sigma / dudvdw... = J * d^n sigma / dxdydz...
// where:
//
// dxdxdz... = J * dudvdw...
//
//                    | dx/du   dx/dv   dx/dw  ... |
//                    |                            |
//     d {x,y,z,...}  | dy/du   dy/dv   dy/dw  ... |
// J = ------------ = |                            |
//     d {u,v,w,...}  | dz/du   dz/dv   dz/dw  ... |
//                    |                            |
//                    |  ...     ...    ...    ... |
//
  SLOG("KineLimits", pDEBUG)
       << "Computing Jacobian for transformation: "
       << KinePhaseSpace::AsString(fromps) << " --> "
       << KinePhaseSpace::AsString(tops);
 
  double J=0;
  bool forward;
  const Kinematics & kine = i->Kine();
 
  // cover the simple case
  if ( fromps == tops )
  {
    forward = true;
    J = 1.;
  }
  //
  // transformation: {Q2}|E -> {lnQ2}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSQ2fE,kPSlogQ2fE,forward) )
  {
    J = 1. / kine.Q2();
  }
 
  //
  // transformation: {QD2}|E -> {Q2}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSQD2fE,kPSQ2fE,forward) )
  {
    J = TMath::Power(1+kine.Q2()/controls::kMQD2,-2)/controls::kMQD2;
  }
 
  //
  // transformation: {x,y}|E -> {lnx,lny}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSxyfE,kPSlogxlogyfE,forward) )
  {
    J = 1. / (kine.x() * kine.y());
  }
 
  //
  // transformation: {log10x,log10Q2}|E -> {x,Q2}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSxQ2fE,kPSlog10xlog10Q2fE,forward) )
  {
    J = TMath::Log(10.)*kine.x() * TMath::Log(10.)*kine.Q2();
  }
 
  //
  // transformation: {Q2,y}|E -> {lnQ2,lny}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSQ2yfE,kPSlogQ2logyfE,forward) )
  {
    J = 1. / (kine.Q2() * kine.y());
  }
 
  //
  // transformation: {x,y}|E -> {x,Q2}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSxQ2fE,kPSxyfE,forward) )
  {
    const InitialState & init_state = i->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M();
    double x  = kine.x();
    J = 2*x*Ev*M;
  }
 
  //
  // transformation: {Q2,y}|E -> {x,y}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSQ2yfE,kPSxyfE,forward) )
  {
    const InitialState & init_state = i->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M();
    double y  = kine.y();
    J = 2*y*Ev*M;
  }
 
  //
  // transformation: {W,Q2}|E -> {W,lnQ2}|E
  //
  else if ( TransformMatched(fromps,tops,kPSWQ2fE,kPSWlogQ2fE,forward) )
  {
    J = 1. / kine.Q2();
  }
 
  //
  // transformation: {W,QD2}|E -> {W,Q2}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSWQD2fE,kPSWQ2fE,forward) )
  {
    J = TMath::Power(1+kine.Q2()/controls::kMQD2,-2)/controls::kMQD2;
  }
 
  //
  // transformation: {W2,Q2}|E --> {x,y}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSW2Q2fE,kPSxyfE,forward) )
  {
    const InitialState & init_state = i->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M();
    double y  = kine.y();
    J = TMath::Power(2*M*Ev,2) * y;
  }
 
  //
  // transformation: {W,Q2}|E -> {x,y}|E
  //
  else
  if ( TransformMatched(fromps,tops,kPSWQ2fE,kPSxyfE,forward) )
  {
    const InitialState & init_state = i->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M();
    double y  = kine.y();
    double W  = kine.W();
    J = 2*TMath::Power(M*Ev,2) * y/W;
  }
 
  // Transformation: {Omegalep,Omegapi}|E -> {Omegalep,Thetapi}|E
  else if ( TransformMatched(fromps,tops,kPSElOlOpifE,kPSElOlTpifE,forward) ) {
    // Use symmetry to turn 4d integral into 3d * 2pi
    J = 2*constants::kPi;
  }
 
  // Transformation: {Tl,ctl} -> {W,Q2}|E
  // IMPORTANT NOTE: This transformation can't be done exactly in two
  // dimensions due to Fermi motion. For fixed {W,Q2}, the azimuthal angle phi
  // is uniformly distributed in the hit nucleon rest frame, while for fixed
  // {Tl,ctl} it is uniform in the lab frame (i.e., the target nucleus rest
  // frame). A 3D Jacobian is needed in order to account for the relationship
  // between these two azimuthal angles. In the implementation below, I choose
  // to neglect Fermi motion. Under this approximation, the lab and hit nucleon
  // rest frames are the same. Doing things this way is a stopgap solution for
  // the new XSecShape_CCMEC reweighting dial developed for MicroBooNE. A full
  // solution should use the Jacobian that connects the 3D phase spaces which
  // each include phi. - S. Gardiner, 29 July 2020
  else if ( TransformMatched(fromps, tops, kPSTlctl, kPSWQ2fE, forward) )
  {
    // Probe properties (mass, energy, momentum)
    const InitialState& init_state = i->InitState();
    double mv = init_state.Probe()->Mass();
    double Ev = init_state.ProbeE( kRfLab );
    double pv = std::sqrt( std::max(0., Ev*Ev - mv*mv) );
 
    // Invariant mass of the initial hit nucleon
    const TLorentzVector& hit_nuc_P4 = init_state.Tgt().HitNucP4();
    double M = hit_nuc_P4.M();
 
    // Outgoing lepton mass
    double ml = i->FSPrimLepton()->Mass();
 
    double W = i->Kine().GetKV( kKVW );
    double Q2 = i->Kine().GetKV( kKVQ2 );
    double El = Ev - ( (W*W + Q2 - M*M) / (2.*M) );
    double pl = std::sqrt( std::max(0., El*El - ml*ml) );
 
    // Compute the Jacobian for {Tl, ctl} --> {W, Q2}
    // (it will be inverted below for the inverse transformation)
    J = W / ( 2. * pv * pl * M );
  }
 
  else {
     std::ostringstream msg;
     msg << "Can not compute Jacobian for transforming: "
       << KinePhaseSpace::AsString(fromps) << " --> "
       << KinePhaseSpace::AsString(tops);
     SLOG("KineLimits", pFATAL) << "*** " << msg.str();
     throw genie::exceptions::InteractionException(msg.str());
     //exit(1);
  }
 
  // if any of the above transforms was reverse, invert the jacobian
  if(!forward) J = 1./J;
 
  return J;
}

References genie::KinePhaseSpace::AsString(), genie::Interaction::FSPrimLepton(), genie::Kinematics::GetKV(), genie::Target::HitNucP4(), genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::Interaction::Kine(), genie::kKVQ2, genie::kKVW, genie::controls::kMQD2, genie::constants::kPi, genie::kPSElOlOpifE, genie::kPSElOlTpifE, genie::kPSlog10xlog10Q2fE, genie::kPSlogQ2fE, genie::kPSlogQ2logyfE, genie::kPSlogxlogyfE, genie::kPSQ2fE, genie::kPSQ2yfE, genie::kPSQD2fE, genie::kPSTlctl, genie::kPSW2Q2fE, genie::kPSWlogQ2fE, genie::kPSWQ2fE, genie::kPSWQD2fE, genie::kPSxQ2fE, genie::kPSxyfE, genie::kRfHitNucRest, genie::kRfLab, pDEBUG, pFATAL, genie::InitialState::Probe(), genie::InitialState::ProbeE(), genie::Kinematics::Q2(), Q2(), SLOG, genie::InitialState::Tgt(), TransformMatched(), genie::Kinematics::W(), W(), genie::Kinematics::x(), and genie::Kinematics::y().

Referenced by genie::RESKinematicsGenerator::ProcessEventRecord(), genie::AhrensDMELPXSec::XSec(), genie::AhrensNCELPXSec::XSec(), genie::AivazisCharmPXSecLO::XSec(), genie::AlvarezRusoCOHPiPXSec::XSec(), genie::BardinIMDRadCorPXSec::XSec(), genie::BergerSehgalCOHPiPXSec2015::XSec(), genie::BostedChristyEMPXSec::XSec(), genie::BSKLNBaseRESPXSec2014::XSec(), genie::DMElectronPXSec::XSec(), genie::EmpiricalMECPXSec2015::XSec(), genie::HEDISPXSec::XSec(), genie::IMDAnnihilationPXSec::XSec(), genie::KovalenkoQELCharmPXSec::XSec(), genie::LwlynSmithQELCCPXSec::XSec(), genie::MKSPPPXSec2020::XSec(), genie::NievesQELCCPXSec::XSec(), genie::NievesSimoVacasMECPXSec2016::XSec(), genie::NuElectronPXSec::XSec(), genie::P33PaschosLalakulichPXSec::XSec(), genie::PaisQELLambdaPXSec::XSec(), genie::QPMDISPXSec::XSec(), genie::QPMDMDISPXSec::XSec(), genie::ReinDFRPXSec::XSec(), genie::ReinSehgalCOHPiPXSec::XSec(), genie::ReinSehgalRESPXSec::XSec(), genie::RosenbluthPXSec::XSec(), genie::SlowRsclCharmDISPXSecLO::XSec(), genie::SmithMonizQELCCPXSec::XSec(), and genie::StrumiaVissaniIBDPXSec::XSec().

PhaseSpaceVolume()

double genie::utils::kinematics::PhaseSpaceVolume	(	const Interaction *const	i,
		KinePhaseSpace_t	ps )

Definition at line 36 of file KineUtils.cxx.

{
  double vol = 0;
 
  const KPhaseSpace & phase_space = in->PhaseSpace();
 
  switch(ps) {
 
  case(kPSQ2fE):
  {
    Range1D_t Q2 = phase_space.Limits(kKVQ2);
    vol = Q2.max - Q2.min;
    return vol;
    break;
  }
  case(kPSq2fE):
  {
    Range1D_t q2 = phase_space.Limits(kKVq2);
    vol = q2.min - q2.max;
    return vol;
    break;
  }
  case(kPSWfE):
  {
    Range1D_t W = phase_space.Limits(kKVW);
    vol = W.max - W.min;
    return vol;
    break;
  }
  case(kPSWQ2fE):
  {
    Range1D_t W = phase_space.Limits(kKVW);
    if(W.max<0) return 0;
    const int kNW = 100;
    double dW = (W.max-W.min)/(kNW-1);
    Interaction interaction(*in);
    for(int iw=0; iw<kNW; iw++) {
      interaction.KinePtr()->SetW(W.min + iw*dW);
      Range1D_t Q2 = interaction.PhaseSpace().Q2Lim_W();
      double dQ2 = (Q2.max-Q2.min);
      vol += (dW*dQ2);
    }
    return vol;
    break;
  }
  case(kPSxyfE):
  {
    Range1D_t W = phase_space.Limits(kKVW);
    if(W.max<0) return 0;
 
    const InitialState & init_state = in->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M();
 
    const int    kNx = 100;
    const int    kNy = 100;
    const double kminx = controls::kMinX;
    const double kminy = controls::kMinY;
    const double kdx = (controls::kMaxX - kminx) / (kNx-1);
    const double kdy = (controls::kMaxY - kminy) / (kNy-1);
    const double kdV = kdx*kdy;
 
    double cW=-1, cQ2 = -1;
 
    Interaction interaction(*in);
 
    for(int ix=0; ix<kNx; ix++) {
      double x = kminx+ix*kdx;
      for(int iy=0; iy<kNy; iy++) {
         double y = kminy+iy*kdy;
 
         XYtoWQ2(Ev, M, cW, cQ2, x, y);
         if(!math::IsWithinLimits(cW, W)) continue;
 
         interaction.KinePtr()->SetW(cW);
         Range1D_t Q2 = interaction.PhaseSpace().Q2Lim_W();
         if(!math::IsWithinLimits(cQ2, Q2)) continue;
 
         vol += kdV;
      }
    }
    return vol;
    break;
  }
  default:
   SLOG("KineLimits", pERROR)
     << "Couldn't compute phase space volume for "
                    << KinePhaseSpace::AsString(ps) << "using \n" << *in;
   return 0;
  }
  return 0;
}

References genie::KinePhaseSpace::AsString(), genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::utils::math::IsWithinLimits(), genie::Interaction::KinePtr(), genie::kKVQ2, genie::kKVq2, genie::kKVW, genie::controls::kMaxX, genie::controls::kMaxY, genie::controls::kMinX, genie::controls::kMinY, genie::kPSQ2fE, genie::kPSq2fE, genie::kPSWfE, genie::kPSWQ2fE, genie::kPSxyfE, genie::kRfHitNucRest, genie::KPhaseSpace::Limits(), genie::Range1D_t::max, genie::Range1D_t::min, pERROR, genie::Interaction::PhaseSpace(), genie::InitialState::ProbeE(), Q2(), genie::KPhaseSpace::Q2Lim_W(), genie::Kinematics::SetW(), SLOG, genie::InitialState::Tgt(), W(), and XYtoWQ2().

Referenced by genie::DFRKinematicsGenerator::ProcessEventRecord(), genie::DISKinematicsGenerator::ProcessEventRecord(), genie::DMDISKinematicsGenerator::ProcessEventRecord(), genie::DMEKinematicsGenerator::ProcessEventRecord(), genie::DMELKinematicsGenerator::ProcessEventRecord(), genie::IBDKinematicsGenerator::ProcessEventRecord(), genie::NuEKinematicsGenerator::ProcessEventRecord(), genie::QELKinematicsGenerator::ProcessEventRecord(), and genie::RESKinematicsGenerator::ProcessEventRecord().

Q2()

double genie::utils::kinematics::Q2

(

const Interaction *const

i

)

Definition at line 1077 of file KineUtils.cxx.

{
// Get Q^2 from kinematics object
// If Q^2 is not set but x,y are, then compute Q2 from x,y
 
  const Kinematics & kinematics = interaction->Kine();
 
  if (kinematics.KVSet(kKVQ2) || kinematics.KVSet(kKVq2)) {
    double Q2 = kinematics.Q2();
    return Q2;
  }
  if (kinematics.KVSet(kKVy)) {
    const InitialState & init_state = interaction->InitState();
    double Mn = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
    double x  = kinematics.x();
    double y  = kinematics.y();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double Q2 = 2*Mn*Ev*x*y;
    return Q2;
  }
  SLOG("KineLimits", pERROR) << "Couldn't compute Q^2 for \n"<< *interaction;
  return 0;
}

References genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::Interaction::Kine(), genie::kKVQ2, genie::kKVq2, genie::kKVy, genie::kRfHitNucRest, pERROR, genie::InitialState::ProbeE(), Q2(), SLOG, and genie::InitialState::Tgt().

Referenced by CEvNSQ2Lim(), CohNuLim(), CohQ2Lim(), Cohq2Lim(), CohW2Lim(), CohYLim(), DarkQ2Lim(), Darkq2Lim(), DarkQ2Lim_W(), Darkq2Lim_W(), genie::utils::kinematics::electromagnetic::InelQ2Lim(), InelQ2Lim(), genie::utils::kinematics::electromagnetic::Inelq2Lim(), Inelq2Lim(), genie::utils::kinematics::electromagnetic::InelQ2Lim_W(), InelQ2Lim_W(), genie::utils::kinematics::electromagnetic::Inelq2Lim_W(), Inelq2Lim_W(), IsAboveCharmThreshold(), Jacobian(), PhaseSpaceVolume(), Q2(), Q2toQD2(), Q2YtoX(), RESImportanceSamplingEnvelope(), SlowRescalingVar(), UpdateWQ2FromXY(), UpdateWYFromXQ2(), UpdateXFromQ2Y(), UpdateXYFromWQ2(), WQ2toXY(), XQ2toWY(), XYtoQ2(), and XYtoWQ2().

Q2toQD2()

double genie::utils::kinematics::Q2toQD2

(

double

Q2

)

Definition at line 1061 of file KineUtils.cxx.

{
// Q2 -> QD2 transformation. QD2 takes out the dipole form of the form factors
// making the differential cross section to be flatter and speeding up the
// kinematical selection.
 
  assert(Q2>0);
  return TMath::Power(1+Q2/controls::kMQD2, -1);
}

References genie::controls::kMQD2, and Q2().

Referenced by genie::RESKinematicsGenerator::ProcessEventRecord().

Q2YtoX()

double genie::utils::kinematics::Q2YtoX	(	double	Ev,
		double	M,
		double	Q2,
		double	y )

Definition at line 1222 of file KineUtils.cxx.

{
// Converts (Q^2,y) => x
// Ev is the neutrino energy at the struck nucleon rest frame
// M is the nucleon mass - it does not need to be on the mass shell
  assert(Ev > 0. && M  > 0. && Q2 > 0. && y  > 0.);
 
  double x = Q2 / (2. * y * M * Ev);
 
  LOG("KineLimits", pDEBUG) << "(Ev=" << Ev << ",Q2=" << Q2
    << ",y=" << y << ",M=" << M << ") => x=" << x;
 
  return x;
}

References LOG, pDEBUG, and Q2().

Referenced by genie::QELEventGeneratorSuSA::SelectLeptonKinematics(), genie::MECGenerator::SelectNSVLeptonKinematics(), genie::MECGenerator::SelectSuSALeptonKinematics(), and UpdateXFromQ2Y().

QD2toQ2()

double genie::utils::kinematics::QD2toQ2

(

double

QD2

)

Definition at line 1071 of file KineUtils.cxx.

{
  assert(QD2>0);
  return controls::kMQD2*(1/QD2-1);
}

References genie::controls::kMQD2.

Referenced by genie::RESKinematicsGenerator::ProcessEventRecord().

RESImportanceSamplingEnvelope()

double genie::utils::kinematics::RESImportanceSamplingEnvelope	(	double *	x,
		double *	par )

Definition at line 1372 of file KineUtils.cxx.

{
  //-- inputs
  double Q2    = x[0];   // Q2
  double W     = x[1];   // invariant mass
 
  //-- parameters
  double mD    = par[0]; // resonance mass
  double gD    = par[1]; // resonance width
  double xsmax = par[2]; // safety factor * max cross section in (W,Q2)
  double Wmax  = par[3]; // kinematically allowed Wmax
  //double E     = par[4]; // neutrino energy
 
//  return 3*xsmax; ////////////////NOTE
 
  xsmax*=5;
 
  double func = 0;
 
  if(Wmax > mD) {
     // -- if the resonance mass is within the kinematical range then
     // -- the envelope is defined as a plateau above the resonance peak,
     // -- a steeply falling leading edge (low-W side) and a more slowly
     // -- falling trailing edge (high-W side)
 
     double hwfe = mD+2*gD; // high W falling edge
     double lwfe = mD-gD/2; // low  W falling edge
 
     if(W < lwfe) {
       //low-W falling edge
       func = xsmax / (1 + 5* TMath::Power((W-lwfe)/gD,2));
     } else if (W > hwfe) {
       //high-W falling edge
       func = xsmax / (1 + TMath::Power((W-hwfe)/gD,2));
     } else {
       // plateau
       func = xsmax;
     }
  } else {
     // -- if the resonance mass is above the kinematical range then the
     // -- envelope is a small plateau just bellow Wmax and a falling edge
     // -- at lower W
 
     double plateau_edge = Wmax-0.1;
 
     if (W > plateau_edge) {
       // plateau
       func = xsmax;
     } else {
       //low-W falling edge
       func = xsmax / (1 + TMath::Power((W-plateau_edge)/gD,2));
     }
  }
 
  if(Q2>controls::kMQD2) {
    func *= (0.5 + 1./(1 + Q2/controls::kMQD2));
  }
 
  return func;
}

References func(), genie::controls::kMQD2, Q2(), and W().

Referenced by genie::RESKinematicsGenerator::LoadConfig().

SlowRescalingVar()

double genie::utils::kinematics::SlowRescalingVar	(	double	x,
		double	Q2,
		double	M,
		double	mc )

Definition at line 1257 of file KineUtils.cxx.

{
// x : scaling variable (plain or modified)
// Q2: momentum transfer
// M : hit nucleon "mass" (nucleon can be off the mass shell)
// mc: charm mass
 
  double mc2 = TMath::Power(mc,2);
  double v   = 0.5*Q2/(M*x);
  double xc  = x + 0.5*mc2/(M*v);
  return xc;
}

References Q2().

Referenced by genie::QPMDISStrucFuncBase::CalcPDFs(), genie::QPMDMDISStrucFuncBase::CalcPDFs(), and IsAboveCharmThreshold().

TransformMatched()

bool genie::utils::kinematics::TransformMatched	(	KinePhaseSpace_t	ia,
		KinePhaseSpace_t	ib,
		KinePhaseSpace_t	a,
		KinePhaseSpace_t	b,
		bool &	fwd )

Definition at line 341 of file KineUtils.cxx.

{
 
  // match a->b or b->a
  bool matched = ( (a==inpa&&b==inpb) || (a==inpb&&b==inpa) );
 
  if(matched) {
    if(a==inpa&&b==inpb) fwd = true;  // matched forward transform: a->b
    else                 fwd = false; // matched reverse transform: b->a
  }
  return matched;
}

References a.

Referenced by Jacobian().

UpdateWQ2FromXY()

void genie::utils::kinematics::UpdateWQ2FromXY

(

const Interaction *

in

)

Definition at line 1290 of file KineUtils.cxx.

{
  Kinematics * kine = in->KinePtr();
 
  if(kine->KVSet(kKVx) && kine->KVSet(kKVy)) {
    const InitialState & init_state = in->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M(); // can be off mass shell
    double x  = kine->x();
    double y  = kine->y();
 
    double Q2=-1,W=-1;
    kinematics::XYtoWQ2(Ev,M,W,Q2,x,y);
    kine->SetQ2(Q2);
    kine->SetW(W);
  }
}

References genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::Interaction::KinePtr(), genie::kKVx, genie::kKVy, genie::kRfHitNucRest, genie::Kinematics::KVSet(), genie::InitialState::ProbeE(), Q2(), genie::Kinematics::SetQ2(), genie::Kinematics::SetW(), genie::InitialState::Tgt(), W(), genie::Kinematics::x(), XYtoWQ2(), and genie::Kinematics::y().

Referenced by genie::DFRKinematicsGenerator::ComputeMaxXSec(), genie::DISKinematicsGenerator::ComputeMaxXSec(), genie::DMDISKinematicsGenerator::ComputeMaxXSec(), genie::utils::gsl::d2XSec_dxdy_E::DoEval(), genie::utils::gsl::d3Xsec_dOmegaldThetapi::DoEval(), genie::utils::gsl::d4Xsec_dEldThetaldOmegapi::DoEval(), genie::utils::gsl::d5Xsec_dEldOmegaldOmegapi::DoEval(), genie::utils::gsl::d5XSecAR::DoEval(), genie::KPhaseSpace::IsAllowed(), main(), genie::DISKinematicsGenerator::ProcessEventRecord(), genie::DMDISKinematicsGenerator::ProcessEventRecord(), genie::COHKinematicsGenerator::SetKinematics(), and genie::KPhaseSpace::TLim().

UpdateWYFromXQ2()

void genie::utils::kinematics::UpdateWYFromXQ2

(

const Interaction *

in

)

Definition at line 1326 of file KineUtils.cxx.

{
  Kinematics * kine = in->KinePtr();
 
  if(kine->KVSet(kKVx) && kine->KVSet(kKVQ2)) {
    const InitialState & init_state = in->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M(); // can be off mass shell
    double x  = kine->x();
    double Q2 = kine->Q2();
 
    double W=-1,y=-1;
    kinematics::XQ2toWY(Ev,M,W,Q2,x,y);
    kine->SetW(W);
    kine->Sety(y);
  }
}

References genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::Interaction::KinePtr(), genie::kKVQ2, genie::kKVx, genie::kRfHitNucRest, genie::Kinematics::KVSet(), genie::InitialState::ProbeE(), genie::Kinematics::Q2(), Q2(), genie::Kinematics::SetW(), genie::Kinematics::Sety(), genie::InitialState::Tgt(), W(), genie::Kinematics::x(), and XQ2toWY().

Referenced by genie::utils::gsl::d2XSec_dlog10xdlog10Q2_E::DoEval(), genie::HEDISKinematicsGenerator::ProcessEventRecord(), and genie::HEDISKinematicsGenerator::Scan().

UpdateXFromQ2Y()

void genie::utils::kinematics::UpdateXFromQ2Y

(

const Interaction *

in

)

Definition at line 1344 of file KineUtils.cxx.

{
  Kinematics * kine = in->KinePtr();
 
  if(kine->KVSet(kKVy) && kine->KVSet(kKVQ2)) {
 
    const ProcessInfo &  pi = in->ProcInfo();
    const InitialState & init_state = in->InitState();
    double M = 0.0;
    double Ev = 0.0;
 
    if (pi.IsCoherentProduction()) {
      M = in->InitState().Tgt().Mass(); // nucleus mass
      Ev = init_state.ProbeE(kRfLab);
    }
    else {
      M = in->InitState().Tgt().HitNucP4Ptr()->M(); //can be off m/shell
      Ev = init_state.ProbeE(kRfHitNucRest);
    }
 
    double y  = kine->y();
    double Q2 = kine->Q2();
 
    double x = kinematics::Q2YtoX(Ev,M,Q2,y);
    kine->Setx(x);
  }
}

References genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherentProduction(), genie::Interaction::KinePtr(), genie::kKVQ2, genie::kKVy, genie::kRfHitNucRest, genie::kRfLab, genie::Kinematics::KVSet(), genie::Target::Mass(), genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::Kinematics::Q2(), Q2(), Q2YtoX(), genie::Kinematics::Setx(), genie::InitialState::Tgt(), and genie::Kinematics::y().

Referenced by genie::COHKinematicsGenerator::CalculateKin_BergerSehgal(), genie::utils::gsl::d2XSec_dQ2dy_E::DoEval(), genie::utils::gsl::d2XSec_dQ2dydt_E::DoEval(), and genie::COHKinematicsGenerator::MaxXSec_BergerSehgal().

UpdateXYFromWQ2()

void genie::utils::kinematics::UpdateXYFromWQ2

(

const Interaction *

in

)

Definition at line 1308 of file KineUtils.cxx.

{
  Kinematics * kine = in->KinePtr();
 
  if(kine->KVSet(kKVW) && kine->KVSet(kKVQ2)) {
    const InitialState & init_state = in->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M(); // can be off mass shell
    double W  = kine->W();
    double Q2 = kine->Q2();
 
    double x=-1,y=-1;
    kinematics::WQ2toXY(Ev,M,W,Q2,x,y);
    kine->Setx(x);
    kine->Sety(y);
  }
}

References genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::Interaction::KinePtr(), genie::kKVQ2, genie::kKVW, genie::kRfHitNucRest, genie::Kinematics::KVSet(), genie::InitialState::ProbeE(), genie::Kinematics::Q2(), Q2(), genie::Kinematics::Setx(), genie::Kinematics::Sety(), genie::InitialState::Tgt(), genie::Kinematics::W(), W(), and WQ2toXY().

W()

double genie::utils::kinematics::W

(

const Interaction *const

i

)

Definition at line 1101 of file KineUtils.cxx.

{
  const ProcessInfo & process_info = interaction->ProcInfo();
 
  if(process_info.IsQuasiElastic()) {
    // hadronic inv. mass is equal to the recoil nucleon on-shell mass
    int rpdgc = interaction->RecoilNucleonPdg();
    double M = PDGLibrary::Instance()->Find(rpdgc)->Mass();
    return M;
  }
 
  const Kinematics & kinematics = interaction->Kine();
  if(kinematics.KVSet(kKVW)) {
    double W = kinematics.W();
    return W;
  }
  if(kinematics.KVSet(kKVx) && kinematics.KVSet(kKVy)) {
    const InitialState & init_state = interaction->InitState();
    double Ev = init_state.ProbeE(kRfHitNucRest);
    double M  = init_state.Tgt().HitNucP4Ptr()->M();
    double M2 = M*M;
    double x  = kinematics.x();
    double y  = kinematics.y();
    double W2 = TMath::Max(0., M2 + 2*Ev*M*y*(1-x));
    double W  = TMath::Sqrt(W2);
    return W;
  }
  SLOG("KineLimits", pERROR) << "Couldn't compute W for \n"<< *interaction;
  return 0;
}

References genie::PDGLibrary::Find(), genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::ProcessInfo::IsQuasiElastic(), genie::Interaction::Kine(), genie::kKVW, genie::kKVx, genie::kKVy, genie::kRfHitNucRest, pERROR, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::Interaction::RecoilNucleonPdg(), SLOG, genie::InitialState::Tgt(), and W().

Referenced by genie::AGKYLowW2019::AssertValidity(), genie::PythiaBaseHadro2019::AssertValidity(), DarkQ2Lim(), DarkQ2Lim_W(), Darkq2Lim_W(), DarkWLim(), genie::AGKYLowW2019::Hadronize(), genie::utils::kinematics::electromagnetic::InelQ2Lim(), InelQ2Lim(), genie::utils::kinematics::electromagnetic::InelQ2Lim_W(), InelQ2Lim_W(), genie::utils::kinematics::electromagnetic::Inelq2Lim_W(), Inelq2Lim_W(), genie::utils::kinematics::electromagnetic::InelWLim(), InelWLim(), IsAboveCharmThreshold(), Jacobian(), genie::AGKYLowW2019::MultiplicityProb(), PhaseSpaceVolume(), genie::KPhaseSpace::Q2Lim_W(), genie::KPhaseSpace::Q2Lim_W_SPP(), genie::KPhaseSpace::Q2Lim_W_SPP_iso(), RESImportanceSamplingEnvelope(), genie::AGKYLowW2019::SelectParticles(), UpdateWQ2FromXY(), UpdateWYFromXQ2(), UpdateXYFromWQ2(), W(), WQ2toXY(), XQ2toWY(), XYtoW(), and XYtoWQ2().

WQ2toXY()

void genie::utils::kinematics::WQ2toXY	(	double	Ev,
		double	M,
		double	W,
		double	Q2,
		double &	x,
		double &	y )

Definition at line 1132 of file KineUtils.cxx.

{
// Converts (W,Q2) => (x,y)
// Uses the system: a) W^2 - M^2 = 2*Ev*M*y*(1-x) and b) Q^2 = 2*x*y*M*Ev
// Ev is the neutrino energy at the struck nucleon rest frame
// M is the nucleon mass - it does not need to be on the mass shell
 
  double M2  = TMath::Power(M,2);
  double W2  = TMath::Power(W,2);
 
  x = Q2 / (W2-M2+Q2);
  y = (W2-M2+Q2) / (2*M*Ev);
 
  x = TMath::Min(1.,x);
  y = TMath::Min(1.,y);
  x = TMath::Max(0.,x);
  y = TMath::Max(0.,y);
 
  LOG("KineLimits", pDEBUG)
        << "(W=" << W << ",Q2=" << Q2 << ") => (x="<< x << ", y=" << y<< ")";
}

References LOG, pDEBUG, Q2(), and W().

Referenced by genie::utils::gsl::d2XSec_dWdQ2_E::DoEval(), genie::DMELEventGenerator::ProcessEventRecord(), genie::DMELKinematicsGenerator::ProcessEventRecord(), genie::IBDKinematicsGenerator::ProcessEventRecord(), genie::QELEventGenerator::ProcessEventRecord(), genie::QELEventGeneratorSM::ProcessEventRecord(), genie::QELKinematicsGenerator::ProcessEventRecord(), genie::RESKinematicsGenerator::ProcessEventRecord(), genie::SPPEventGenerator::ProcessEventRecord(), genie::MECGenerator::SelectEmpiricalKinematics(), genie::DMELKinematicsGenerator::SpectralFuncExperimentalCode(), genie::QELKinematicsGenerator::SpectralFuncExperimentalCode(), UpdateXYFromWQ2(), and genie::EmpiricalMECPXSec2015::XSec().

XQ2toWY()

void genie::utils::kinematics::XQ2toWY	(	double	Ev,
		double	M,
		double &	W,
		double	Q2,
		double	x,
		double &	y )

Definition at line 1173 of file KineUtils.cxx.

{
// Converts (x,Q2) => (W,Y)
// Uses the system: a) W^2 - M^2 = 2*Ev*M*y*(1-x) and b) Q^2 = 2*x*y*M*Ev
// Ev is the neutrino energy at the struck nucleon rest frame
// M is the nucleon mass - it does not need to be on the mass shell
 
  double M2  = TMath::Power(M,2);
  y = Q2 / (2 * x * M * Ev);
  y = TMath::Min(1.,y);
 
  double W2  = M2 + 2*Ev*M*y*(1-x);
  W  = TMath::Sqrt(TMath::Max(0., W2));
 
  LOG("KineLimits", pDEBUG)
      << "(x=" << x << ",Q2=" << Q2 << " => (W=" << W << ",Y=" << y << ")";
}

References LOG, pDEBUG, Q2(), and W().

Referenced by UpdateWYFromXQ2().

XYtoQ2()

double genie::utils::kinematics::XYtoQ2	(	double	Ev,
		double	M,
		double	x,
		double	y )

Definition at line 1208 of file KineUtils.cxx.

{
// Converts (x,y) => Q2
// Ev is the neutrino energy at the struck nucleon rest frame
// M is the nucleon mass - it does not need to be on the mass shell
 
  double Q2 = 2*x*y*M*Ev;
 
  LOG("KineLimits", pDEBUG) << "(x=" << x << ",y=" << y << ") => Q2=" << Q2;
 
  return Q2;
}

References LOG, pDEBUG, and Q2().

Referenced by genie::QPMDISPXSec::XSec(), and genie::QPMDMDISPXSec::XSec().

XYtoW()

double genie::utils::kinematics::XYtoW	(	double	Ev,
		double	M,
		double	x,
		double	y )

Definition at line 1192 of file KineUtils.cxx.

{
// Converts (x,y) => W
// Ev is the neutrino energy at the struck nucleon rest frame
// M is the nucleon mass - it does not need to be on the mass shell
 
  double M2  = TMath::Power(M,2);
  double W2  = M2 + 2*Ev*M*y*(1-x);
  double W  = TMath::Sqrt(TMath::Max(0., W2));
 
  LOG("KineLimits", pDEBUG) << "(x=" << x << ",y=" << y << ") => W=" << W;
 
  return W;
}

References LOG, pDEBUG, and W().

Referenced by genie::KNOTunedQPMDISPXSec::DISRESJoinSuppressionFactor(), genie::QELEventGeneratorSuSA::SelectLeptonKinematics(), genie::MECGenerator::SelectNSVLeptonKinematics(), and genie::MECGenerator::SelectSuSALeptonKinematics().

XYtoWQ2()

void genie::utils::kinematics::XYtoWQ2	(	double	Ev,
		double	M,
		double &	W,
		double &	Q2,
		double	x,
		double	y )

Definition at line 1155 of file KineUtils.cxx.

{
// Converts (x,y) => (W,Q2)
// Uses the system: a) W^2 - M^2 = 2*Ev*M*y*(1-x) and b) Q^2 = 2*x*y*M*Ev
// Ev is the neutrino energy at the struck nucleon rest frame
// M is the nucleon mass - it does not need to be on the mass shell
 
  double M2  = TMath::Power(M,2);
  double W2  = M2 + 2*Ev*M*y*(1-x);
 
  W  = TMath::Sqrt(TMath::Max(0., W2));
  Q2 = 2*x*y*M*Ev;
 
  LOG("KineLimits", pDEBUG)
      << "(x=" << x << ",y=" << y << " => (W=" << W << ",Q2=" << Q2 << ")";
}

References LOG, pDEBUG, Q2(), and W().

Referenced by GetCrossSection(), PhaseSpaceVolume(), genie::DFRKinematicsGenerator::ProcessEventRecord(), genie::DISKinematicsGenerator::ProcessEventRecord(), genie::DMDISKinematicsGenerator::ProcessEventRecord(), and UpdateWQ2FromXY().