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 mpi, 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 mpi, 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 mpi, double mlep, double Ev)
Range1D_t	Cohq2Lim (double Mn, double mpi, double mlep, 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 mpi)
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)
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	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 <costas.andreopoulos stfc.ac.uk> University of Liverpool & STFC Rutherford Appleton Lab

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

November 26, 2004

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

Function Documentation

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

Definition at line 1166 of file KineUtils.cxx.

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

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

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

Definition at line 1343 of file KineUtils.cxx.

References d, and func().

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

{
  //-- 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 745 of file KineUtils.cxx.

References genie::Range1D_t::max, genie::Range1D_t::min, and ana::Sqrt().

Referenced by CohYLim().

{
  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	mpi,
		double	mlep,
		double	Ev
	)

Definition at line 673 of file KineUtils.cxx.

References genie::units::A, a, b, C, plot_validation_datamc::c, CohW2Min(), cet::sqlite::max(), genie::Range1D_t::max, genie::Range1D_t::min, pERROR, Q2(), SLOG, and ana::Sqrt().

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

{
  // The expressions for Q^2 min appears in PRD 74, 054007 (2006) by
  // Kartavtsev, Paschos, and Gounaris

  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, mpi);

  // 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	mpi,
		double	mlep,
		double	Ev
	)

Definition at line 711 of file KineUtils.cxx.

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

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

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

Definition at line 720 of file KineUtils.cxx.

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

Referenced by CohYLim().

{
  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
  // Kartavtsev, Paschos, and Gounaris

  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, mpi);
  W2l.max = (T1 - T2 * T2 ) *
    (1.0 / Mnterm) *
    (1.0 / (Q2 + mlep * mlep));

  return W2l;
}

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

Definition at line 798 of file KineUtils.cxx.

Referenced by CohQ2Lim(), and CohW2Lim().

{
  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
  // Kartavtsev, Paschos, and Gounaris

  return (Mn + mpi) * (Mn + mpi);
}

Range1D_t genie::utils::kinematics::CohXLim

(

void

)

Definition at line 665 of file KineUtils.cxx.

References genie::controls::kASmallNum, and submit_syst::x.

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

{
// 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	mpi,
		double	mlep,
		double	Ev,
		double	Q2,
		double	xsi
	)

Definition at line 762 of file KineUtils.cxx.

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

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

{
  Range1D_t ylim;
  ylim.min = -1;
  ylim.max = -1;

  Range1D_t W2lim = genie::utils::kinematics::CohW2Lim(Mn, mpi, 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 << "; mpi = " << mpi << "; 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 787 of file KineUtils.cxx.

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

{
// 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 887 of file KineUtils.cxx.

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

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

{
// 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 903 of file KineUtils.cxx.

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

{
// 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 829 of file KineUtils.cxx.

References ana::assert(), genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, Q2(), SLOG, ana::Sqrt(), and W.

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

{
// 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 875 of file KineUtils.cxx.

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

{
// 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 806 of file KineUtils.cxx.

References ana::assert(), genie::controls::kASmallNum, genie::constants::kNeutronMass, genie::constants::kPhotontest, genie::Range1D_t::max, genie::Range1D_t::min, ana::Sqrt(), and W.

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

{
// 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 915 of file KineUtils.cxx.

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

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

{
// 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 935 of file KineUtils.cxx.

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

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

{
  // 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 951 of file KineUtils.cxx.

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

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

{
  // 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 1311 of file KineUtils.cxx.

References func(), xmax, and make_mec_shifts_plots::xmin.

{
  //-- 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 359 of file KineUtils.cxx.

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

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

{
// 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 375 of file KineUtils.cxx.

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

{
// 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 309 of file KineUtils.cxx.

References ana::assert(), genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, Q2(), SLOG, and ana::Sqrt().

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

{
// 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 347 of file KineUtils.cxx.

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

{
// 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 288 of file KineUtils.cxx.

References ana::assert(), genie::controls::kASmallNum, genie::constants::kNeutronMass, genie::constants::kPhotontest, genie::Range1D_t::max, genie::Range1D_t::min, ana::Sqrt(), and W.

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

{
// 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 387 of file KineUtils.cxx.

References ana::assert(), genie::controls::kASmallNum, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, SLOG, and submit_syst::x.

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

{
// 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::kASmallNum;
  x.max = 1.         - controls::kASmallNum;

  return x;
}

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

Definition at line 407 of file KineUtils.cxx.

References ana::assert(), MECModelEnuComparisons::i, InelXLim(), InelYLim_X(), genie::controls::kASmallNum, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, SLOG, submit_syst::x, and submit_syst::y.

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

{
// 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::kASmallNum);
    y.max = TMath::Min(y.max, 1 - controls::kASmallNum);
  } 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 442 of file KineUtils.cxx.

References genie::units::A, a, ana::assert(), b, plot_validation_datamc::c, d, genie::controls::kASmallNum, genie::Range1D_t::max, genie::Range1D_t::min, pDEBUG, SLOG, ana::Sqrt(), submit_syst::x, and submit_syst::y.

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

{
// 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::kASmallNum;
  y.max = TMath::Min(1., A+B) - controls::kASmallNum;

  return y;
}

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

Definition at line 1133 of file KineUtils.cxx.

References genie::constants::kLightestChmHad, SlowRescalingVar(), ana::Sqrt(), registry_explorer::v, W, Wmin, and submit_syst::x.

Referenced by genie::QPMDISStrucFuncBase::CalcPDFs().

{
// 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 128 of file KineUtils.cxx.

References genie::KinePhaseSpace::AsString(), exit(), genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::utils::mec::J(), genie::Interaction::Kine(), genie::controls::kMQD2, genie::constants::kPi, genie::kPSElOlOpifE, genie::kPSElOlTpifE, genie::kPSlogQ2fE, genie::kPSlogQ2logyfE, genie::kPSlogxlogyfE, genie::kPSQ2fE, genie::kPSQ2yfE, genie::kPSQD2fE, genie::kPSW2Q2fE, genie::kPSWlogQ2fE, genie::kPSWQ2fE, genie::kPSWQD2fE, genie::kPSxyfE, genie::kRfHitNucRest, pDEBUG, pFATAL, genie::InitialState::ProbeE(), genie::Kinematics::Q2(), SLOG, genie::InitialState::Tgt(), TransformMatched(), W, genie::Kinematics::W(), genie::Kinematics::x(), genie::Kinematics::y(), and submit_syst::y.

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

{
// 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();

  //
  // transformation: {Q2}|E -> {lnQ2}|E
  //
  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: {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;
  }

  else {
     SLOG("KineLimits", pFATAL)
       << "*** Can not compute Jacobian for transforming: "
       << KinePhaseSpace::AsString(fromps) << " --> "
       << KinePhaseSpace::AsString(tops);
     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 34 of file KineUtils.cxx.

References genie::KinePhaseSpace::AsString(), genie::Target::HitNucP4Ptr(), in, 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(), q2, genie::KPhaseSpace::Q2Lim_W(), genie::Kinematics::SetW(), SLOG, genie::InitialState::Tgt(), W, submit_syst::x, XYtoWQ2(), and submit_syst::y.

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

{
  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 991 of file KineUtils.cxx.

References genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::Interaction::Kine(), kinematics(), genie::kKVQ2, genie::kKVq2, genie::kKVy, genie::kRfHitNucRest, genie::Kinematics::KVSet(), pERROR, genie::InitialState::ProbeE(), genie::Kinematics::Q2(), SLOG, genie::InitialState::Tgt(), genie::Kinematics::x(), submit_syst::x, genie::Kinematics::y(), and submit_syst::y.

Referenced by genie::BaryonResonanceDecayer::AcceptPionDecay(), genie::MECGenerator::AddFinalStateLepton(), genie::GRV98LO::AllPDFs(), genie::BBA03ELFormFactorsModel::BBA03Fit(), BuildStdNtuple(), ana::CalcMECDoubleGaussEnhShiftedParam(), ana::CalcMinosResSupp(), CalculateFormFactor(), genie::COHHadronicSystemGenerator::CalculateHadronicSystem_BergerSehgalFM(), genie::COHKinematicsGenerator::CalculateKin_AlvarezRuso(), genie::SKKinematicsGenerator::CalculateKin_AtharSingleKaon(), genie::COHKinematicsGenerator::CalculateKin_BergerSehgal(), genie::COHKinematicsGenerator::CalculateKin_BergerSehgalFM(), CohQ2Lim(), Cohq2Lim(), CohYLim(), genie::utils::ComputeFullQELPXSec(), genie::QELEventGeneratorSM::ComputeMaxDiffv(), genie::RESKinematicsGenerator::ComputeMaxXSec(), genie::IBDKinematicsGenerator::ComputeMaxXSec(), genie::QELKinematicsGenerator::ComputeMaxXSec(), genie::DMELKinematicsGenerator::ComputeMaxXSec(), genie::QELEventGeneratorSM::ComputeMaxXSec(), genie::QELEventGeneratorSM::ComputeMaxXSec2(), novarwgt::ConvertGenieEvent(), ConvertToGST(), DarkQ2Lim(), Darkq2Lim(), DarkQ2Lim_W(), Darkq2Lim_W(), genie::utils::gsl::dXSec_dQ2_E::DoEval(), genie::utils::gsl::d2Xsec_dQ2dv::DoEval(), genie::utils::gsl::d2XSecRESFast_dWQ2_E::DoEval(), genie::utils::gsl::d2XSec_dQ2dy_E::DoEval(), genie::utils::gsl::d2XSec_dQ2dydt_E::DoEval(), genie::utils::gsl::d2XSec_dWdQ2_E::DoEval(), genie::utils::gsl::d2XSec_dWdQ2_EW::DoEval(), genie::utils::gsl::d5XSecAR::DoEval(), genie::utils::gsl::d5Xsec_dEldOmegaldOmegapi::DoEval(), genie::utils::gsl::d4Xsec_dEldThetaldOmegapi::DoEval(), genie::utils::gsl::d3Xsec_dOmegaldThetapi::DoEval(), genie::KovalenkoQELCharmPXSec::DR(), genie::SmithMonizQELCCPXSec::dsQES_dQ2_SM(), genie::BergerSehgalFMCOHPiPXSec2015::ExactKinematicTerm(), genie::BergerSehgalCOHPiPXSec2015::ExactKinematicTerm(), genie::LwlynSmithQELCCPXSec::FullDifferentialXSec(), GetCrossSection(), genie::TransverseEnhancementFFModel::GetTransEnhMagFF(), GetTunedMECWeight(), genie::LHAPDF6::Gluon(), InelQ2Lim(), Inelq2Lim(), genie::utils::kinematics::electromagnetic::InelQ2Lim(), genie::utils::kinematics::electromagnetic::Inelq2Lim(), InelQ2Lim_W(), Inelq2Lim_W(), genie::utils::kinematics::electromagnetic::InelQ2Lim_W(), genie::utils::kinematics::electromagnetic::Inelq2Lim_W(), genie::GRV98LO::Initialize(), ana::NuWROSyst::InitializeHistograms(), genie::KPhaseSpace::IsAllowed(), ana::kMECDoubleGaussEnh(), main(), supernova::SnovaGen::MakeMCTruth(), MakePlots(), genie::COHKinematicsGenerator::MaxXSec_BergerSehgal(), genie::COHKinematicsGenerator::MaxXSec_BergerSehgalFM(), genie::COHKinematicsGenerator::MaxXSec_ReinSehgal(), PhaseSpaceVolume(), genie::SmithMonizUtils::PhaseSpaceVolume(), genie::BergerSehgalFMCOHPiPXSec2015::PionCOMAbsMomentum(), genie::BergerSehgalCOHPiPXSec2015::PionCOMAbsMomentum(), genie::NuEPrimaryLeptonGenerator::ProcessEventRecord(), genie::PrimaryLeptonGenerator::ProcessEventRecord(), genie::COHElKinematicsGenerator::ProcessEventRecord(), genie::QELKinematicsGenerator::ProcessEventRecord(), genie::OutgoingDarkGenerator::ProcessEventRecord(), genie::RESKinematicsGenerator::ProcessEventRecord(), genie::IBDKinematicsGenerator::ProcessEventRecord(), genie::DMELKinematicsGenerator::ProcessEventRecord(), genie::KPhaseSpace::q2Lim(), genie::KPhaseSpace::q2Lim_W(), genie::utils::phys::R99118(), SaveToRootFile(), genie::MECGenerator::SelectEmpiricalKinematics(), genie::MECGenerator::SelectNSVLeptonKinematics(), genie::COHKinematicsGenerator::SetKinematics(), ana::MECDoubleGaussEnhSyst::Shift(), ana::MinosResSuppSyst::Shift(), ana::MECDoubleGaussEnhSystUP::Shift(), ana::MECDoubleGaussEnhSystDOWN::Shift(), genie::QELKinematicsGenerator::SpectralFuncExperimentalCode(), genie::DMELKinematicsGenerator::SpectralFuncExperimentalCode(), genie::KPhaseSpace::TLim(), UpdateWQ2FromXY(), UpdateXFromQ2Y(), UpdateXYFromWQ2(), validation(), WQ2toXY(), genie::KovalenkoQELCharmPXSec::xiBar(), genie::COHElasticPXSec::XSec(), genie::EmpiricalMECPXSec2015::XSec(), genie::SlowRsclCharmDISPXSecLO::XSec(), genie::AhrensNCELPXSec::XSec(), genie::ReinDFRPXSec::XSec(), genie::RosenbluthPXSec::XSec(), genie::QPMDISPXSec::XSec(), genie::AivazisCharmPXSecLO::XSec(), genie::AhrensDMELPXSec::XSec(), genie::P33PaschosLalakulichPXSec::XSec(), genie::BergerSehgalFMCOHPiPXSec2015::XSec(), genie::QPMDMDISPXSec::XSec(), genie::BergerSehgalCOHPiPXSec2015::XSec(), genie::ReinSehgalCOHPiPXSec::XSec(), genie::KovalenkoQELCharmPXSec::XSec(), genie::ReinSehgalRESPXSec::XSec(), genie::BSKLNBaseRESPXSec2014::XSec(), genie::NievesQELCCPXSec::XSec(), XYtoQ2(), and genie::KPhaseSpace::YLim().

{
// 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 975 of file KineUtils.cxx.

References ana::assert(), and genie::controls::kMQD2.

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

{
// 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 1117 of file KineUtils.cxx.

References ana::assert(), LOG, pDEBUG, and submit_syst::x.

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

{
// 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 985 of file KineUtils.cxx.

References ana::assert(), and genie::controls::kMQD2.

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

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

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

Definition at line 1249 of file KineUtils.cxx.

References func(), W, and Wmax.

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

{
  //-- 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 1152 of file KineUtils.cxx.

References registry_explorer::v, and submit_syst::x.

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

{
// 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 271 of file KineUtils.cxx.

Referenced by Jacobian().

{

  // 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 1185 of file KineUtils.cxx.

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(), submit_syst::x, XYtoWQ2(), genie::Kinematics::y(), and submit_syst::y.

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

{
  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::UpdateXFromQ2Y

(

const Interaction *

in

)

Definition at line 1221 of file KineUtils.cxx.

References genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::ProcessInfo::IsCoherent(), genie::Interaction::KinePtr(), genie::kKVQ2, genie::kKVy, genie::kRfHitNucRest, genie::kRfLab, genie::Kinematics::KVSet(), genie::Target::Mass(), python.hepunit::pi, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::Kinematics::Q2(), Q2(), Q2YtoX(), genie::Kinematics::Setx(), genie::InitialState::Tgt(), submit_syst::x, genie::Kinematics::y(), and submit_syst::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().

{
  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.IsCoherent()) {
      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 1203 of file KineUtils.cxx.

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(), W, genie::Kinematics::W(), WQ2toXY(), submit_syst::x, and submit_syst::y.

{
  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 1015 of file KineUtils.cxx.

References genie::PDGLibrary::Find(), genie::Target::HitNucP4Ptr(), genie::Interaction::InitState(), genie::PDGLibrary::Instance(), genie::ProcessInfo::IsQuasiElastic(), genie::Interaction::Kine(), kinematics(), genie::kKVW, genie::kKVx, genie::kKVy, genie::kRfHitNucRest, genie::Kinematics::KVSet(), pERROR, genie::InitialState::ProbeE(), genie::Interaction::ProcInfo(), genie::Interaction::RecoilNucleonPdg(), SLOG, ana::Sqrt(), genie::InitialState::Tgt(), W, genie::Kinematics::W(), genie::Kinematics::x(), submit_syst::x, genie::Kinematics::y(), and submit_syst::y.

Referenced by genie::PythiaHadronization::AssertValidity(), genie::KNOHadronization::AssertValidity(), genie::KNOHadronization::Hadronize(), genie::KNOHadronization::MultiplicityProb(), and genie::KNOHadronization::SelectParticles().

{
  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 1046 of file KineUtils.cxx.

References LOG, pDEBUG, and Q2().

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

{
// 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<< ")";
}

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

Definition at line 1103 of file KineUtils.cxx.

References LOG, pDEBUG, and Q2().

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

{
// 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 1087 of file KineUtils.cxx.

References LOG, pDEBUG, ana::Sqrt(), W, and submit_syst::x.

Referenced by genie::QPMDISPXSec::DISRESJoinSuppressionFactor(), genie::QPMDMDISPXSec::DMDISRESJoinSuppressionFactor(), and genie::MECGenerator::SelectNSVLeptonKinematics().

{
// 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 1069 of file KineUtils.cxx.

References LOG, pDEBUG, ana::Sqrt(), and submit_syst::x.

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

{
// 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 << ")";
}