genie::utils::kinematics Namespace Reference

Kinematical utilities. More...

Namespaces
	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 <constantinos.andreopoulos cern.ch> University of Liverpool & STFC Rutherford Appleton Laboratory

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

November 26, 2004

Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org

Function Documentation

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

Definition at line 1258 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;
  }
}

Range1D_t genie::utils::kinematics::CEvNSQ2Lim

(

double

Ev

)

Definition at line 873 of file KineUtils.cxx.

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

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

Definition at line 1453 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;
}

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

Definition at line 810 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;
}

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

Definition at line 730 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;
}

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

Definition at line 772 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;
}

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

Definition at line 781 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;
}

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

Definition at line 861 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);
}

Range1D_t genie::utils::kinematics::CohXLim

(

void

)

Definition at line 722 of file KineUtils.cxx.

{
// Computes x limits for coherent v interactions

  Range1D_t x(controls::kASmallNum, 1.-controls::kASmallNum);
  return x;
}

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

Definition at line 827 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;
}

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

Definition at line 852 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;
}

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

Definition at line 960 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;
}

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

Definition at line 976 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;
}

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

Definition at line 902 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;
}

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

Definition at line 948 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;
}

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

Definition at line 879 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;
}

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

Definition at line 988 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;
}

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

Definition at line 1008 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;
}

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

Definition at line 1024 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;
}

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

Definition at line 1421 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;
}

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

Definition at line 416 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;
}

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

Definition at line 432 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;
}

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

Definition at line 366 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;
}

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

Definition at line 404 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;
}

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

Definition at line 345 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;
}

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

Definition at line 444 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;
}

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

Definition at line 464 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;
}

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

Definition at line 499 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;
}

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

Definition at line 1225 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;
}

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: {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;
}

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;
}

double genie::utils::kinematics::Q2

(

const Interaction *const

i

)

Definition at line 1064 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;
}

double genie::utils::kinematics::Q2toQD2

(

double

Q2

)

Definition at line 1048 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);
}

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

Definition at line 1209 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;
}

double genie::utils::kinematics::QD2toQ2

(

double

QD2

)

Definition at line 1058 of file KineUtils.cxx.

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

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

Definition at line 1359 of file KineUtils.cxx.

{
  //-- inputs
  double QD2   = x[0];   // QD2 (Q2 transformed to take out the dipole form)
  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(QD2<0.5) {
    func *= (1 - (0.5-QD2));
  }

  return func;
}

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

Definition at line 1244 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;
}

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

Definition at line 328 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;
}

void genie::utils::kinematics::UpdateWQ2FromXY

(

const Interaction *

in

)

Definition at line 1277 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);
  }
}

void genie::utils::kinematics::UpdateWYFromXQ2

(

const Interaction *

in

)

Definition at line 1313 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);
  }
}

void genie::utils::kinematics::UpdateXFromQ2Y

(

const Interaction *

in

)

Definition at line 1331 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);
  }
}

void genie::utils::kinematics::UpdateXYFromWQ2

(

const Interaction *

in

)

Definition at line 1295 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);
  }
}

double genie::utils::kinematics::W

(

const Interaction *const

i

)

Definition at line 1088 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;
}

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

Definition at line 1119 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<< ")";
}

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

Definition at line 1160 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 << ")";
}

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

Definition at line 1195 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;
}

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

Definition at line 1179 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;
}

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

Definition at line 1142 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 << ")";
}