Namespaces | Functions
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.

1168 {
1169  // if the min,max are within the existing limits, the cut can be applied
1170  // by narrowing down the xisting limits
1171  if ( utils::math::IsWithinLimits(min_cut, range ) ) range.min = min_cut;
1172  if ( utils::math::IsWithinLimits(max_cut, range ) ) range.max = max_cut;
1173 
1174  // if the min-cut is above the existing max-limit or
1175  // if the max-cut is below the existing min-limit then
1176  // the range should be invalidated
1177 
1178  if (min_cut > range.max || max_cut < range.min) {
1179 
1180  range.min = 0;
1181  range.max = 0;
1182  }
1183 }
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:265
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().

1345 {
1346  //-- inputs
1347  double xb = x[0]; // x
1348  double yb = x[1]; // y
1349 
1350  //-- parameters
1351  double xsmax = 3*par[0]; // safety factor * max cross section in (x,y)
1352  double Ev = par[1]; // neutrino energy;
1353 
1354  if(yb<0.|| yb>1.) return 0.;
1355  if(xb<0.|| xb>1.) return 0.;
1356 
1357  if(Ev<1) return xsmax;
1358  if(xb/Ev<1E-4 && yb>0.95) return 5*xsmax;
1359 
1360  double func = 0;
1361  double xp = 0.1;
1362  double yp = (Ev>2.5) ? 2.5/Ev : 1;
1363 
1364  if(xb>xp) {
1365  double xs0=0;
1366  if(yb<yp) {
1367  xs0 = xsmax;
1368  } else {
1369  xs0 = xsmax - (yb-yp)*xsmax;
1370  }
1371  double d = TMath::Power( (xb-xp)/0.075, 2);
1372  func = xs0/(1 + d);
1373  } else {
1374  if(yb>yp) {
1375  func = xsmax - (yb-yp)*xsmax;
1376  } else {
1377  func = xsmax;
1378  }
1379  }
1380  return func;
1381 }
Int_t par
Definition: SimpleIterate.C:24
Float_t d
Definition: plot.C:236
double func(double x, double y)
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().

747 {
748  Range1D_t nul;
749  nul.min = -1;
750  nul.max = -1;
751 
752  double nu_min = (W2min + Q2 - Mn * Mn) / (2.0 * Mn);
753  double nu_max = (W2max + Q2 - Mn * Mn) / (2.0 * Mn);
754  double xsiQ = xsi * TMath::Sqrt(Q2);
755 
756  nul.min = (xsiQ > nu_min) ? xsiQ : nu_min;
757  nul.max = nu_max;
758 
759  return nul;
760 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double min
Definition: Range1.h:53
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().

674 {
675  // The expressions for Q^2 min appears in PRD 74, 054007 (2006) by
676  // Kartavtsev, Paschos, and Gounaris
677 
678  Range1D_t Q2;
679  Q2.min = 0.0;
680  Q2.max = std::numeric_limits<double>::max(); // Value must be overriden in user options
681 
682  double Mn2 = Mn * Mn;
683  double mlep2 = mlep * mlep;
684  double s = Mn2 + 2.0 * Mn * Ev;
685  double W2min = CohW2Min(Mn, mpi);
686 
687  // Looks like Q2min = A * B - C, where A, B, and C are complicated
688  double a = 1.0;
689  double b = mlep2 / s;
690  double c = W2min / s;
691  double lambda = a * a + b * b + c * c - 2.0 * a * b - 2.0 * a * c - 2.0 * b * c;
692  if (lambda > 0) {
693  double A = (s - Mn * Mn) / 2.0;
694  double B = 1 - TMath::Sqrt(lambda);
695  double C = 0.5 * (W2min + mlep2 - Mn2 * (W2min - mlep2) / s );
696  if (A * B - C < 0) {
697  SLOG("KineLimits", pERROR)
698  << "Q2 kinematic limits calculation failed for CohQ2Lim. "
699  << "Assuming Q2min = 0.0";
700  }
701  Q2.min = TMath::Max(0., A * B - C);
702  } else {
703  SLOG("KineLimits", pERROR)
704  << "Q2 kinematic limits calculation failed for CohQ2Lim. "
705  << "Assuming Q2min = 0.0";
706  }
707 
708  return Q2;
709 }
#define pERROR
Definition: Messenger.h:60
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double CohW2Min(double Mn, double mpi)
Definition: KineUtils.cxx:798
const double C
const XML_Char * s
Definition: expat.h:262
const double a
static const double A
Definition: Units.h:82
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const hit & b
Definition: hits.cxx:21
double min
Definition: Range1.h:53
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
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.

712 {
713  Range1D_t Q2 = utils::kinematics::CohQ2Lim(Mn, mpi, mlep, Ev);
714  Range1D_t q2;
715  q2.min = - Q2.max;
716  q2.max = - Q2.min;
717  return q2;
718 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Double_t q2[12][num]
Definition: f2_nu.C:137
Range1D_t CohQ2Lim(double Mn, double mpi, double mlep, double Ev)
Definition: KineUtils.cxx:673
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
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().

722 {
723  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
724  // Kartavtsev, Paschos, and Gounaris
725 
726  Range1D_t W2l;
727  W2l.min = -1;
728  W2l.max = -1;
729 
730  double s = Mn * Mn + 2.0 * Mn * Ev;
731  double Mnterm = 1 - Mn * Mn / s;
732  double Mlterm = 1 - mlep * mlep / s;
733  // Here T1, T2 are generically "term 1" and "term 2" in a long expression
734  double T1 = 0.25 * s * s * Mnterm * Mnterm * Mlterm;
735  double T2 = Q2 - (0.5 * s * Mnterm) + (0.5 * mlep * mlep * Mnterm);
736 
737  W2l.min = CohW2Min(Mn, mpi);
738  W2l.max = (T1 - T2 * T2 ) *
739  (1.0 / Mnterm) *
740  (1.0 / (Q2 + mlep * mlep));
741 
742  return W2l;
743 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double CohW2Min(double Mn, double mpi)
Definition: KineUtils.cxx:798
const XML_Char * s
Definition: expat.h:262
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
double genie::utils::kinematics::CohW2Min ( double  Mn,
double  mpi 
)

Definition at line 798 of file KineUtils.cxx.

Referenced by CohQ2Lim(), and CohW2Lim().

799 {
800  // These expressions for W^2 min and max appear in PRD 74, 054007 (2006) by
801  // Kartavtsev, Paschos, and Gounaris
802 
803  return (Mn + mpi) * (Mn + mpi);
804 }
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().

666 {
667 // Computes x limits for coherent v interactions
668 
670  return x;
671 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
static const double kASmallNum
Definition: Controls.h:40
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().

764 {
765  Range1D_t ylim;
766  ylim.min = -1;
767  ylim.max = -1;
768 
769  Range1D_t W2lim = genie::utils::kinematics::CohW2Lim(Mn, mpi, mlep, Ev, Q2);
770  if (W2lim.min > W2lim.max) {
771  LOG("KineLimits", pDEBUG)
772  << "Kinematically forbidden region in CohYLim. W2min = " << W2lim.min
773  << "; W2max =" << W2lim.max;
774  LOG("KineLimits", pDEBUG)
775  << " Mn = " << Mn << "; mpi = " << mpi << "; mlep = "
776  << mlep << "; Ev = " << Ev << "; Q2 = " << Q2;
777  return ylim;
778  }
780  Q2, Mn, xsi);
781  ylim.min = nulim.min / Ev;
782  ylim.max = nulim.max / Ev;
783 
784  return ylim;
785 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Range1D_t CohNuLim(double W2min, double W2max, double Q2, double Mn, double xsi)
Definition: KineUtils.cxx:745
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
Range1D_t CohW2Lim(double Mn, double mpi, double mlep, double Ev, double Q2)
Definition: KineUtils.cxx:720
#define pDEBUG
Definition: Messenger.h:64
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.

788 {
789 // Computes y limits for coherent v interactions
790 
792  1.-ml/EvL - controls::kASmallNum);
793  return y;
794 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
static const double kASmallNum
Definition: Controls.h:40
static const double kPionMass
Definition: Constants.h:74
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().

889 {
890 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
891 
892  Range1D_t Q2;
893  Q2.min = -1;
894  Q2.max = -1;
895 
897  if(W.min<0) return Q2;
898 
899  Q2 = utils::kinematics::DarkQ2Lim_W(Ev,M,ml,W.min,Q2min_cut);
900  return Q2;
901 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:829
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
#define W(x)
Range1D_t DarkWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:806
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.

905 {
906 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
907 
908  Range1D_t Q2 = utils::kinematics::DarkQ2Lim(Ev,M,ml,-1.*q2min_cut);
909  Range1D_t q2;
910  q2.min = - Q2.max;
911  q2.max = - Q2.min;
912  return q2;
913 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Double_t q2[12][num]
Definition: f2_nu.C:137
Range1D_t DarkQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:887
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
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().

831 {
832 // Computes Q2 limits (>0) @ the input W for inelastic v interactions
833 
834  Range1D_t Q2;
835  Q2.min = -1;
836  Q2.max = -1;
837 
838  double M2 = TMath::Power(M, 2.);
839  double ml2 = TMath::Power(ml, 2.);
840  double W2 = TMath::Power(W, 2.);
841  double s = M2 + 2*M*Ev + ml2;
842  assert(s > 0);
843  double sqs = TMath::Sqrt(s);
844  double E1CM = (s + ml2 - M2) / (2.*sqs);
845  double p1CM = TMath::Max(0., E1CM*E1CM - ml2);
846  p1CM = TMath::Sqrt(p1CM);
847  double E3CM = (s + ml2 - W2) / (2.*sqs);
848  double p3CM = TMath::Max(0., E3CM*E3CM - ml2);
849  p3CM = TMath::Sqrt(p3CM);
850 
851  SLOG("KineLimits", pDEBUG) << "s = " << s;
852  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
853  SLOG("KineLimits", pDEBUG) << "M = " << M;
854  SLOG("KineLimits", pDEBUG) << "W = " << W;
855  SLOG("KineLimits", pDEBUG) << "E1_CM = " << E1CM;
856  SLOG("KineLimits", pDEBUG) << "p1_CM = " << p1CM;
857  SLOG("KineLimits", pDEBUG) << "E3_CM = " << E3CM;
858  SLOG("KineLimits", pDEBUG) << "p3_CM = " << p3CM;
859 
860  Q2.min = TMath::Power(p3CM - p1CM,2) - TMath::Power((W2 - M2) / (2.*sqs),2);
861  Q2.max = TMath::Power(p3CM + p1CM,2) - TMath::Power((W2 - M2) / (2.*sqs),2);
862 
863  SLOG("KineLimits", pDEBUG) << "Nominal Q^2 limits: " << Q2.min << " , " << Q2.max;
864  // guard against overflows
865  Q2.max = TMath::Max(0., Q2.max);
866  Q2.min = TMath::Max(0., Q2.min);
867 
868  // limit the minimum Q2
869  if(Q2.min < Q2min_cut) {Q2.min = Q2min_cut; }
870  if(Q2.max < Q2.min ) {Q2.min = -1; Q2.max = -1;}
871 
872  return Q2;
873 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
const XML_Char * s
Definition: expat.h:262
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
#define W(x)
#define pDEBUG
Definition: Messenger.h:64
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.

877 {
878 // Computes q2 (<0) limits @ the input W for inelastic v interactions
879 
880  Range1D_t Q2 = utils::kinematics::DarkQ2Lim_W(Ev,M,ml,W,-1.*q2min_cut);
881  Range1D_t q2;
882  q2.min = - Q2.max;
883  q2.max = - Q2.min;
884  return q2;
885 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:829
Double_t q2[12][num]
Definition: f2_nu.C:137
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
#define W(x)
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().

807 {
808 // Computes W limits for inelastic v interactions
809 //
810  double M2 = TMath::Power(M,2);
811  double ml2 = TMath::Power(ml,2);
812  double s = M2 + 2*M*Ev + ml2;
813  assert (s>0);
814 
815  Range1D_t W;
817 // W.min = kNeutronMass + kPionMass;
818  W.max = TMath::Sqrt(s) - ml;
819  if(W.max<=W.min) {
820  W.min = -1;
821  W.max = -1;
822  return W;
823  }
826  return W;
827 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
static const double kPhotontest
Definition: Constants.h:80
const XML_Char * s
Definition: expat.h:262
static const double kNeutronMass
Definition: Constants.h:77
static const double kASmallNum
Definition: Controls.h:40
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
#define W(x)
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().

917 {
918 // Computes Bjorken x limits for inelastic interactions
919 // For the dark matter case, it is relatively straightforward
920  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
921  double Wmin = Wl.min;
922  double W2min = Wmin*Wmin;
923  SLOG("KineLimits", pDEBUG) << "W^2_min = " << W2min;
924  Range1D_t Q2l = utils::kinematics::DarkQ2Lim_W(Ev, M, ml, Wmin);
925  SLOG("KineLimits", pDEBUG) << "Q^2 range : " << Q2l.min << " , " << Q2l.max;
926  double M2 = M*M;
927  Range1D_t x;
928  x.min = Q2l.min / (Q2l.min + W2min - M2);
929  x.max = Q2l.max / (Q2l.max + W2min - M2);
930 
931  SLOG("KineLimits", pDEBUG) << "x = [" << x.min << ", " << x.max << "]";
932  return x;
933 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:829
double max
Definition: Range1.h:54
double Wmin[kNWBins]
double min
Definition: Range1.h:53
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
Range1D_t DarkWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:806
#define pDEBUG
Definition: Messenger.h:64
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().

936 {
937  // For dark inelastic scattering, can compute exactly and is much simpler
938  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
939  double Wmin = Wl.min;
940  double W2min = Wmin*Wmin;
941  Range1D_t Q2l = utils::kinematics::DarkQ2Lim_W(Ev, M, ml, Wmin);
942  double M2 = M*M;
943  Range1D_t y;
944  y.min = (Q2l.min + W2min - M2) / (2*Ev*M);
945  y.max = (Q2l.max + W2min - M2) / (2*Ev*M);
946 
947  SLOG("KineLimits", pDEBUG) << "y = [" << y.min << ", " << y.max << "]";
948  return y;
949 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:829
double max
Definition: Range1.h:54
double Wmin[kNWBins]
double min
Definition: Range1.h:53
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
Range1D_t DarkWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:806
#define pDEBUG
Definition: Messenger.h:64
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().

954 {
955  // Computes y limits @ the input x for inelastic interactions
956  // We hit y_min when W = W_min and y_max when Q^2 = Q_min^2 or Q^2 = Q_max^2
957 
958  Range1D_t y;
959  y.min = -1;
960  y.max = -1;
961 
962  Range1D_t Wl = utils::kinematics::DarkWLim(Ev, M, ml);
963  double Wmin = Wl.min;
964  double W2min = Wmin*Wmin;
965  double M2 = M*M;
966  double ml2 = ml*ml;
967  y.min = (W2min - M2) / (1.-x) / (2*Ev*M);
968  y.max = 2.* M * x *(Ev*Ev - ml2) / Ev / (2. * M * Ev * x + M2 * x * x + ml2);
969 
970  return y;
971 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
double max
Definition: Range1.h:54
double Wmin[kNWBins]
double min
Definition: Range1.h:53
Range1D_t DarkWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:806
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.

1313 {
1314  //-- inputs
1315  double xb = x[0]; // scaling variable x brorken
1316  //double y = x[1]; // inelasticity y
1317 
1318  //-- parameters
1319  double xpeak = par[0]; // x peak
1320  double xmin = par[1]; // min x
1321  double xmax = 1.; // max x
1322  double xsmax = par[2]; // safety factor * max cross section in (x,y)
1323 
1324  double func = 0;
1325 
1326  if(xb < xpeak/20.) {
1327  //low-x falling edge
1328  double plateau_edge = xpeak/20.;
1329  double slope = xsmax/(xmin - plateau_edge);
1330  func = xsmax - slope * (xb - plateau_edge);
1331  } else if (xb > 2*xpeak) {
1332  //high-x falling edge
1333  double plateau_edge = 2*xpeak;
1334  double slope = xsmax/(xmax - plateau_edge);
1335  func = xsmax - slope * (xb - plateau_edge);
1336  } else {
1337  // plateau
1338  func = xsmax;
1339  }
1340  return func;
1341 }
std::map< std::string, double > xmax
Int_t par
Definition: SimpleIterate.C:24
double func(double x, double y)
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().

361 {
362 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
363 
364  Range1D_t Q2;
365  Q2.min = -1;
366  Q2.max = -1;
367 
369  if(W.min<0) return Q2;
370 
371  Q2 = utils::kinematics::InelQ2Lim_W(Ev,M,ml,W.min,Q2min_cut);
372  return Q2;
373 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
Range1D_t InelWLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:288
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:309
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
#define W(x)
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.

377 {
378 // Computes Q2 (>0) limits irrespective of W for inelastic v interactions
379 
380  Range1D_t Q2 = utils::kinematics::InelQ2Lim(Ev,M,ml,-1.*q2min_cut);
381  Range1D_t q2;
382  q2.min = - Q2.max;
383  q2.max = - Q2.min;
384  return q2;
385 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t InelQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:359
Double_t q2[12][num]
Definition: f2_nu.C:137
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
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().

311 {
312 // Computes Q2 limits (>0) @ the input W for inelastic v interactions
313 
314  Range1D_t Q2;
315  Q2.min = -1;
316  Q2.max = -1;
317 
318  double M2 = TMath::Power(M, 2.);
319  double ml2 = TMath::Power(ml, 2.);
320  double W2 = TMath::Power(W, 2.);
321  double s = M2 + 2*M*Ev;
322 
323  SLOG("KineLimits", pDEBUG) << "s = " << s;
324  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
325  assert (s>0);
326 
327  double auxC = 0.5*(s-M2)/s;
328  double aux1 = s + ml2 - W2;
329  double aux2 = aux1*aux1 - 4*s*ml2;
330 
331  (aux2 < 0) ? ( aux2 = 0 ) : ( aux2 = TMath::Sqrt(aux2) );
332 
333  Q2.max = -ml2 + auxC * (aux1 + aux2); // => 0
334  Q2.min = -ml2 + auxC * (aux1 - aux2); // => 0
335 
336  // guard against overflows
337  Q2.max = TMath::Max(0., Q2.max);
338  Q2.min = TMath::Max(0., Q2.min);
339 
340  // limit the minimum Q2
341  if(Q2.min < Q2min_cut) {Q2.min = Q2min_cut; }
342  if(Q2.max < Q2.min ) {Q2.min = -1; Q2.max = -1;}
343 
344  return Q2;
345 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
const XML_Char * s
Definition: expat.h:262
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
#define W(x)
#define pDEBUG
Definition: Messenger.h:64
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.

349 {
350 // Computes q2 (<0) limits @ the input W for inelastic v interactions
351 
352  Range1D_t Q2 = utils::kinematics::InelQ2Lim_W(Ev,M,ml,W,-1.*q2min_cut);
353  Range1D_t q2;
354  q2.min = - Q2.max;
355  q2.max = - Q2.min;
356  return q2;
357 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
Definition: KineUtils.cxx:309
Double_t q2[12][num]
Definition: f2_nu.C:137
double max
Definition: Range1.h:54
double min
Definition: Range1.h:53
#define W(x)
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().

289 {
290 // Computes W limits for inelastic v interactions
291 //
292  double M2 = TMath::Power(M,2);
293  double s = M2 + 2*M*Ev;
294  assert (s>0);
295 
296  Range1D_t W;
298  W.max = TMath::Sqrt(s) - ml;
299  if(W.max<=W.min) {
300  W.min = -1;
301  W.max = -1;
302  return W;
303  }
306  return W;
307 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
static const double kPhotontest
Definition: Constants.h:80
const XML_Char * s
Definition: expat.h:262
static const double kNeutronMass
Definition: Constants.h:77
static const double kASmallNum
Definition: Controls.h:40
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
#define W(x)
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().

389 {
390 // Computes Bjorken x limits for inelastic v interactions
391 
392  double M2 = TMath::Power(M, 2.);
393  double ml2 = TMath::Power(ml,2.);
394  double s = M2 + 2*M*Ev;
395 
396  SLOG("KineLimits", pDEBUG) << "s = " << s;
397  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
398  assert (s>M2);
399 
400  Range1D_t x;
401  x.min = ml2/(s-M2) + controls::kASmallNum;
402  x.max = 1. - controls::kASmallNum;
403 
404  return x;
405 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
const XML_Char * s
Definition: expat.h:262
static const double kASmallNum
Definition: Controls.h:40
double max
Definition: Range1.h:54
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
#define pDEBUG
Definition: Messenger.h:64
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().

408 {
409 // Computes y limits for inelastic v interactions
410 
411  Range1D_t y;
412  y.min = 999;
413  y.max = -999;
414 
415  Range1D_t xl = kinematics::InelXLim(Ev,M,ml);
416  assert(xl.min>0 && xl.max>0);
417 
418  const unsigned int N=100;
419  const double logxmin = TMath::Log10(xl.min);
420  const double logxmax = TMath::Log10(xl.max);
421  const double dlogx = (logxmax-logxmin) / (double)(N-1);
422 
423  for(unsigned int i=0; i<N; i++) {
424  double x = TMath::Power(10, logxmin + i*dlogx);
425 
426  Range1D_t y_x = kinematics::InelYLim_X(Ev,M,ml,x);
427  if(y_x.min>=0 && y_x.min<=1) y.min = TMath::Min(y.min, y_x.min);
428  if(y_x.max>=0 && y_x.max<=1) y.max = TMath::Max(y.max, y_x.max);
429  }
430 
431  if(y.max >= 0 && y.max <= 1 && y.min >= 0 && y.min <= 1) {
432  y.min = TMath::Max(y.min, controls::kASmallNum);
433  y.max = TMath::Min(y.max, 1 - controls::kASmallNum);
434  } else {
435  y.min = -1;
436  y.max = -1;
437  }
438  SLOG("KineLimits", pDEBUG) << "y = [" << y.min << ", " << y.max << "]";
439  return y;
440 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
Range1D_t InelXLim(double Ev, double M, double ml)
Definition: KineUtils.cxx:387
static const double kASmallNum
Definition: Controls.h:40
double max
Definition: Range1.h:54
Range1D_t InelYLim_X(double Ev, double M, double ml, double x)
Definition: KineUtils.cxx:442
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
#define pDEBUG
Definition: Messenger.h:64
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().

445 {
446 // Computes y limits @ the input x for inelastic v interactions
447 
448  Range1D_t y;
449  y.min = -1;
450  y.max = -1;
451 
452  double Ev2 = TMath::Power(Ev,2);
453  double ml2 = TMath::Power(ml,2);
454 
455  SLOG("KineLimits", pDEBUG) << "x = " << x;
456  SLOG("KineLimits", pDEBUG) << "Ev = " << Ev;
457 
458  assert (Ev>0);
459  assert (x>0&&x<1);
460 
461  double a = 0.5 * ml2/(M*Ev*x);
462  double b = ml2/Ev2;
463  double c = 1 + 0.5*x*M/Ev;
464  double d = TMath::Max(0., TMath::Power(1-a,2.) - b);
465 
466  double A = 0.5 * (1-a-0.5*b)/c;
467  double B = 0.5 * TMath::Sqrt(d)/c;
468 
469  y.min = TMath::Max(0., A-B) + controls::kASmallNum;
470  y.max = TMath::Min(1., A+B) - controls::kASmallNum;
471 
472  return y;
473 }
A simple [min,max] interval for doubles.
Definition: Range1.h:43
const double a
Float_t d
Definition: plot.C:236
static const double kASmallNum
Definition: Controls.h:40
static const double A
Definition: Units.h:82
double max
Definition: Range1.h:54
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const hit & b
Definition: hits.cxx:21
assert(nhit_max >=nhit_nbins)
double min
Definition: Range1.h:53
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
#define pDEBUG
Definition: Messenger.h:64
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().

1135 {
1136 // x : scaling variable (plain or modified)
1137 // Q2: momentum transfer
1138 // M : hit nucleon "mass" (nucleon can be off the mass shell)
1139 // mc: charm mass
1140 //
1141  double M2 = TMath::Power(M,2);
1142  double v = 0.5*Q2/(M*x);
1143  double W2 = TMath::Max(0., M2+2*M*v-Q2);
1144  double W = TMath::Sqrt(W2);
1145  double Wmin = M + kLightestChmHad;
1146  double xc = utils::kinematics::SlowRescalingVar(x,Q2,M,mc);
1147 
1148  if(xc>=1 || W<=Wmin) return false;
1149  else return true;
1150 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
static const double kLightestChmHad
Definition: Constants.h:79
double SlowRescalingVar(double x, double Q2, double M, double mc)
Definition: KineUtils.cxx:1152
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
double Wmin[kNWBins]
#define W(x)
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().

130 {
131 // Returns the Jacobian for a kinematical transformation:
132 // from_ps (x,y,z,...) -> to_ps (u,v,w,...).
133 //
134 // Note on the convention used here:
135 // The differential cross-section d^n sigma / dxdydz..., is transformed to
136 // d^n sigma / dudvdw... using the Jacobian, computed in this function, as:
137 // as d^n sigma / dudvdw... = J * d^n sigma / dxdydz...
138 // where:
139 //
140 // dxdxdz... = J * dudvdw...
141 //
142 // | dx/du dx/dv dx/dw ... |
143 // | |
144 // d {x,y,z,...} | dy/du dy/dv dy/dw ... |
145 // J = ------------ = | |
146 // d {u,v,w,...} | dz/du dz/dv dz/dw ... |
147 // | |
148 // | ... ... ... ... |
149 //
150  SLOG("KineLimits", pDEBUG)
151  << "Computing Jacobian for transformation: "
152  << KinePhaseSpace::AsString(fromps) << " --> "
153  << KinePhaseSpace::AsString(tops);
154 
155  double J=0;
156  bool forward;
157  const Kinematics & kine = i->Kine();
158 
159  //
160  // transformation: {Q2}|E -> {lnQ2}|E
161  //
162  if ( TransformMatched(fromps,tops,kPSQ2fE,kPSlogQ2fE,forward) )
163  {
164  J = 1. / kine.Q2();
165  }
166 
167  //
168  // transformation: {QD2}|E -> {Q2}|E
169  //
170  else
171  if ( TransformMatched(fromps,tops,kPSQD2fE,kPSQ2fE,forward) )
172  {
173  J = TMath::Power(1+kine.Q2()/controls::kMQD2,-2)/controls::kMQD2;
174  }
175 
176  //
177  // transformation: {x,y}|E -> {lnx,lny}|E
178  //
179  else
180  if ( TransformMatched(fromps,tops,kPSxyfE,kPSlogxlogyfE,forward) )
181  {
182  J = 1. / (kine.x() * kine.y());
183  }
184 
185  //
186  // transformation: {Q2,y}|E -> {lnQ2,lny}|E
187  //
188  else
189  if ( TransformMatched(fromps,tops,kPSQ2yfE,kPSlogQ2logyfE,forward) )
190  {
191  J = 1. / (kine.Q2() * kine.y());
192  }
193 
194  //
195  // transformation: {Q2,y}|E -> {x,y}|E
196  //
197  else
198  if ( TransformMatched(fromps,tops,kPSQ2yfE,kPSxyfE,forward) )
199  {
200  const InitialState & init_state = i->InitState();
201  double Ev = init_state.ProbeE(kRfHitNucRest);
202  double M = init_state.Tgt().HitNucP4Ptr()->M();
203  double y = kine.y();
204  J = 2*y*Ev*M;
205  }
206 
207  //
208  // transformation: {W,Q2}|E -> {W,lnQ2}|E
209  //
210  else if ( TransformMatched(fromps,tops,kPSWQ2fE,kPSWlogQ2fE,forward) )
211  {
212  J = 1. / kine.Q2();
213  }
214 
215  //
216  // transformation: {W,QD2}|E -> {W,Q2}|E
217  //
218  else
219  if ( TransformMatched(fromps,tops,kPSWQD2fE,kPSWQ2fE,forward) )
220  {
221  J = TMath::Power(1+kine.Q2()/controls::kMQD2,-2)/controls::kMQD2;
222  }
223 
224  //
225  // transformation: {W2,Q2}|E --> {x,y}|E
226  //
227  else
228  if ( TransformMatched(fromps,tops,kPSW2Q2fE,kPSxyfE,forward) )
229  {
230  const InitialState & init_state = i->InitState();
231  double Ev = init_state.ProbeE(kRfHitNucRest);
232  double M = init_state.Tgt().HitNucP4Ptr()->M();
233  double y = kine.y();
234  J = TMath::Power(2*M*Ev,2) * y;
235  }
236 
237  //
238  // transformation: {W,Q2}|E -> {x,y}|E
239  //
240  else
241  if ( TransformMatched(fromps,tops,kPSWQ2fE,kPSxyfE,forward) )
242  {
243  const InitialState & init_state = i->InitState();
244  double Ev = init_state.ProbeE(kRfHitNucRest);
245  double M = init_state.Tgt().HitNucP4Ptr()->M();
246  double y = kine.y();
247  double W = kine.W();
248  J = 2*TMath::Power(M*Ev,2) * y/W;
249  }
250 
251  // Transformation: {Omegalep,Omegapi}|E -> {Omegalep,Thetapi}|E
252  else if ( TransformMatched(fromps,tops,kPSElOlOpifE,kPSElOlTpifE,forward) ) {
253  // Use symmetry to turn 4d integral into 3d * 2pi
254  J = 2*constants::kPi;
255  }
256 
257  else {
258  SLOG("KineLimits", pFATAL)
259  << "*** Can not compute Jacobian for transforming: "
260  << KinePhaseSpace::AsString(fromps) << " --> "
261  << KinePhaseSpace::AsString(tops);
262  exit(1);
263  }
264 
265  // if any of the above transforms was reverse, invert the jacobian
266  if(!forward) J = 1./J;
267 
268  return J;
269 }
const double kPi
double W(bool selected=false) const
Definition: Kinematics.cxx:167
bool TransformMatched(KinePhaseSpace_t ia, KinePhaseSpace_t ib, KinePhaseSpace_t a, KinePhaseSpace_t b, bool &fwd)
Definition: KineUtils.cxx:271
double J(double q0, double q3, double Enu, double ml)
Definition: MECUtils.cxx:141
#define pFATAL
Definition: Messenger.h:57
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double x(bool selected=false) const
Definition: Kinematics.cxx:109
double y(bool selected=false) const
Definition: Kinematics.cxx:122
const Kinematics & Kine(void) const
Definition: Interaction.h:71
static const double kMQD2
Definition: Controls.h:65
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
exit(0)
const InitialState & InitState(void) const
Definition: Interaction.h:69
const char * AsString(Resonance_t res)
resonance id -> string
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
const Target & Tgt(void) const
Definition: InitialState.h:67
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
#define W(x)
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
#define pDEBUG
Definition: Messenger.h:64
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().

36 {
37  double vol = 0;
38 
39  const KPhaseSpace & phase_space = in->PhaseSpace();
40 
41  switch(ps) {
42 
43  case(kPSQ2fE):
44  {
45  Range1D_t Q2 = phase_space.Limits(kKVQ2);
46  vol = Q2.max - Q2.min;
47  return vol;
48  break;
49  }
50  case(kPSq2fE):
51  {
52  Range1D_t q2 = phase_space.Limits(kKVq2);
53  vol = q2.min - q2.max;
54  return vol;
55  break;
56  }
57  case(kPSWfE):
58  {
59  Range1D_t W = phase_space.Limits(kKVW);
60  vol = W.max - W.min;
61  return vol;
62  break;
63  }
64  case(kPSWQ2fE):
65  {
66  Range1D_t W = phase_space.Limits(kKVW);
67  if(W.max<0) return 0;
68  const int kNW = 100;
69  double dW = (W.max-W.min)/(kNW-1);
70  Interaction interaction(*in);
71  for(int iw=0; iw<kNW; iw++) {
72  interaction.KinePtr()->SetW(W.min + iw*dW);
73  Range1D_t Q2 = interaction.PhaseSpace().Q2Lim_W();
74  double dQ2 = (Q2.max-Q2.min);
75  vol += (dW*dQ2);
76  }
77  return vol;
78  break;
79  }
80  case(kPSxyfE):
81  {
82  Range1D_t W = phase_space.Limits(kKVW);
83  if(W.max<0) return 0;
84 
85  const InitialState & init_state = in->InitState();
86  double Ev = init_state.ProbeE(kRfHitNucRest);
87  double M = init_state.Tgt().HitNucP4Ptr()->M();
88 
89  const int kNx = 100;
90  const int kNy = 100;
91  const double kminx = controls::kMinX;
92  const double kminy = controls::kMinY;
93  const double kdx = (controls::kMaxX - kminx) / (kNx-1);
94  const double kdy = (controls::kMaxY - kminy) / (kNy-1);
95  const double kdV = kdx*kdy;
96 
97  double cW=-1, cQ2 = -1;
98 
99  Interaction interaction(*in);
100 
101  for(int ix=0; ix<kNx; ix++) {
102  double x = kminx+ix*kdx;
103  for(int iy=0; iy<kNy; iy++) {
104  double y = kminy+iy*kdy;
105 
106  XYtoWQ2(Ev, M, cW, cQ2, x, y);
107  if(!math::IsWithinLimits(cW, W)) continue;
108 
109  interaction.KinePtr()->SetW(cW);
110  Range1D_t Q2 = interaction.PhaseSpace().Q2Lim_W();
111  if(!math::IsWithinLimits(cQ2, Q2)) continue;
112 
113  vol += kdV;
114  }
115  }
116  return vol;
117  break;
118  }
119  default:
120  SLOG("KineLimits", pERROR)
121  << "Couldn't compute phase space volume for "
122  << KinePhaseSpace::AsString(ps) << "using \n" << *in;
123  return 0;
124  }
125  return 0;
126 }
static const double kMaxX
Definition: Controls.h:44
#define pERROR
Definition: Messenger.h:60
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
A simple [min,max] interval for doubles.
Definition: Range1.h:43
bool IsWithinLimits(double x, Range1D_t range)
Definition: MathUtils.cxx:265
static const double kMinY
Definition: Controls.h:45
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
Double_t q2[12][num]
Definition: f2_nu.C:137
Summary information for an interaction.
Definition: Interaction.h:56
const Var kMaxY([](const caf::SRProxy *sr){if(!sr->vtx.elastic.IsValid) return-1000.0f;if(sr->vtx.elastic.fuzzyk.nshwlid==0) return-1000.0f;float maxyall=-999.0;for(unsigned int i=0;i< sr->vtx.elastic.fuzzyk.nshwlid;i++){maxyall=std::max(std::max(sr->vtx.elastic.fuzzyk.png[i].shwlid.start.Y(), sr->vtx.elastic.fuzzyk.png[i].shwlid.stop.Y()), maxyall);}return maxyall;})
Definition: NueVars.h:60
void XYtoWQ2(double Ev, double M, double &W, double &Q2, double x, double y)
Definition: KineUtils.cxx:1069
Kinematical phase space.
Definition: KPhaseSpace.h:34
static const double kMinX
Definition: Controls.h:43
double max
Definition: Range1.h:54
ifstream in
Definition: comparison.C:7
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
const char * AsString(Resonance_t res)
resonance id -> string
double min
Definition: Range1.h:53
const Target & Tgt(void) const
Definition: InitialState.h:67
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
#define W(x)
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
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(), 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::QELKinematicsGenerator::ProcessEventRecord(), genie::COHElKinematicsGenerator::ProcessEventRecord(), genie::RESKinematicsGenerator::ProcessEventRecord(), genie::OutgoingDarkGenerator::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::MECDoubleGaussEnhSystUP::Shift(), ana::MECDoubleGaussEnhSystDOWN::Shift(), ana::MinosResSuppSyst::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::AivazisCharmPXSecLO::XSec(), genie::QPMDISPXSec::XSec(), genie::P33PaschosLalakulichPXSec::XSec(), genie::AhrensDMELPXSec::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().

992 {
993 // Get Q^2 from kinematics object
994 // If Q^2 is not set but x,y are, then compute Q2 from x,y
995 
996  const Kinematics & kinematics = interaction->Kine();
997 
998  if (kinematics.KVSet(kKVQ2) || kinematics.KVSet(kKVq2)) {
999  double Q2 = kinematics.Q2();
1000  return Q2;
1001  }
1002  if (kinematics.KVSet(kKVy)) {
1003  const InitialState & init_state = interaction->InitState();
1004  double Mn = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell
1005  double x = kinematics.x();
1006  double y = kinematics.y();
1007  double Ev = init_state.ProbeE(kRfHitNucRest);
1008  double Q2 = 2*Mn*Ev*x*y;
1009  return Q2;
1010  }
1011  SLOG("KineLimits", pERROR) << "Couldn't compute Q^2 for \n"<< *interaction;
1012  return 0;
1013 }
#define pERROR
Definition: Messenger.h:60
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double x(bool selected=false) const
Definition: Kinematics.cxx:109
double y(bool selected=false) const
Definition: Kinematics.cxx:122
bool KVSet(KineVar_t kv) const
Definition: Kinematics.cxx:327
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
const Target & Tgt(void) const
Definition: InitialState.h:67
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
void kinematics()
Definition: kinematics.C:10
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
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().

976 {
977 // Q2 -> QD2 transformation. QD2 takes out the dipole form of the form factors
978 // making the differential cross section to be flatter and speeding up the
979 // kinematical selection.
980 
981  assert(Q2>0);
982  return TMath::Power(1+Q2/controls::kMQD2, -1);
983 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
static const double kMQD2
Definition: Controls.h:65
assert(nhit_max >=nhit_nbins)
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().

1119 {
1120 // Converts (Q^2,y) => x
1121 // Ev is the neutrino energy at the struck nucleon rest frame
1122 // M is the nucleon mass - it does not need to be on the mass shell
1123  assert(Ev > 0. && M > 0. && Q2 > 0. && y > 0.);
1124 
1125  double x = Q2 / (2. * y * M * Ev);
1126 
1127  LOG("KineLimits", pDEBUG) << "(Ev=" << Ev << ",Q2=" << Q2
1128  << ",y=" << y << ",M=" << M << ") => x=" << x;
1129 
1130  return x;
1131 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
assert(nhit_max >=nhit_nbins)
#define pDEBUG
Definition: Messenger.h:64
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().

986 {
987  assert(QD2>0);
988  return controls::kMQD2*(1/QD2-1);
989 }
static const double kMQD2
Definition: Controls.h:65
assert(nhit_max >=nhit_nbins)
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().

1251 {
1252  //-- inputs
1253  double QD2 = x[0]; // QD2 (Q2 transformed to take out the dipole form)
1254  double W = x[1]; // invariant mass
1255 
1256  //-- parameters
1257  double mD = par[0]; // resonance mass
1258  double gD = par[1]; // resonance width
1259  double xsmax = par[2]; // safety factor * max cross section in (W,Q2)
1260  double Wmax = par[3]; // kinematically allowed Wmax
1261  //double E = par[4]; // neutrino energy
1262 
1263 // return 3*xsmax; ////////////////NOTE
1264 
1265  xsmax*=5;
1266 
1267  double func = 0;
1268 
1269  if(Wmax > mD) {
1270  // -- if the resonance mass is within the kinematical range then
1271  // -- the envelope is defined as a plateau above the resonance peak,
1272  // -- a steeply falling leading edge (low-W side) and a more slowly
1273  // -- falling trailing edge (high-W side)
1274 
1275  double hwfe = mD+2*gD; // high W falling edge
1276  double lwfe = mD-gD/2; // low W falling edge
1277 
1278  if(W < lwfe) {
1279  //low-W falling edge
1280  func = xsmax / (1 + 5* TMath::Power((W-lwfe)/gD,2));
1281  } else if (W > hwfe) {
1282  //high-W falling edge
1283  func = xsmax / (1 + TMath::Power((W-hwfe)/gD,2));
1284  } else {
1285  // plateau
1286  func = xsmax;
1287  }
1288  } else {
1289  // -- if the resonance mass is above the kinematical range then the
1290  // -- envelope is a small plateau just bellow Wmax and a falling edge
1291  // -- at lower W
1292 
1293  double plateau_edge = Wmax-0.1;
1294 
1295  if (W > plateau_edge) {
1296  // plateau
1297  func = xsmax;
1298  } else {
1299  //low-W falling edge
1300  func = xsmax / (1 + TMath::Power((W-plateau_edge)/gD,2));
1301  }
1302  }
1303 
1304  if(QD2<0.5) {
1305  func *= (1 - (0.5-QD2));
1306  }
1307 
1308  return func;
1309 }
Int_t par
Definition: SimpleIterate.C:24
double Wmax[kNWBins]
double func(double x, double y)
#define W(x)
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().

1154 {
1155 // x : scaling variable (plain or modified)
1156 // Q2: momentum transfer
1157 // M : hit nucleon "mass" (nucleon can be off the mass shell)
1158 // mc: charm mass
1159 
1160  double mc2 = TMath::Power(mc,2);
1161  double v = 0.5*Q2/(M*x);
1162  double xc = x + 0.5*mc2/(M*v);
1163  return xc;
1164 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
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().

274 {
275 
276  // match a->b or b->a
277  bool matched = ( (a==inpa&&b==inpb) || (a==inpb&&b==inpa) );
278 
279  if(matched) {
280  if(a==inpa&&b==inpb) fwd = true; // matched forward transform: a->b
281  else fwd = false; // matched reverse transform: b->a
282  }
283  return matched;
284 }
const double a
const hit & b
Definition: hits.cxx:21
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().

1186 {
1187  Kinematics * kine = in->KinePtr();
1188 
1189  if(kine->KVSet(kKVx) && kine->KVSet(kKVy)) {
1190  const InitialState & init_state = in->InitState();
1191  double Ev = init_state.ProbeE(kRfHitNucRest);
1192  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off mass shell
1193  double x = kine->x();
1194  double y = kine->y();
1195 
1196  double Q2=-1,W=-1;
1197  kinematics::XYtoWQ2(Ev,M,W,Q2,x,y);
1198  kine->SetQ2(Q2);
1199  kine->SetW(W);
1200  }
1201 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
void SetQ2(double Q2, bool selected=false)
Definition: Kinematics.cxx:265
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double x(bool selected=false) const
Definition: Kinematics.cxx:109
double y(bool selected=false) const
Definition: Kinematics.cxx:122
bool KVSet(KineVar_t kv) const
Definition: Kinematics.cxx:327
void XYtoWQ2(double Ev, double M, double &W, double &Q2, double x, double y)
Definition: KineUtils.cxx:1069
void SetW(double W, bool selected=false)
Definition: Kinematics.cxx:289
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
const InitialState & InitState(void) const
Definition: Interaction.h:69
const Target & Tgt(void) const
Definition: InitialState.h:67
#define W(x)
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
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().

1222 {
1223  Kinematics * kine = in->KinePtr();
1224 
1225  if(kine->KVSet(kKVy) && kine->KVSet(kKVQ2)) {
1226 
1227  const ProcessInfo & pi = in->ProcInfo();
1228  const InitialState & init_state = in->InitState();
1229  double M = 0.0;
1230  double Ev = 0.0;
1231 
1232  if (pi.IsCoherent()) {
1233  M = in->InitState().Tgt().Mass(); // nucleus mass
1234  Ev = init_state.ProbeE(kRfLab);
1235  }
1236  else {
1237  M = in->InitState().Tgt().HitNucP4Ptr()->M(); //can be off m/shell
1238  Ev = init_state.ProbeE(kRfHitNucRest);
1239  }
1240 
1241  double y = kine->y();
1242  double Q2 = kine->Q2();
1243 
1244  double x = kinematics::Q2YtoX(Ev,M,Q2,y);
1245  kine->Setx(x);
1246  }
1247 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double Mass(void) const
Definition: Target.cxx:241
double y(bool selected=false) const
Definition: Kinematics.cxx:122
bool KVSet(KineVar_t kv) const
Definition: Kinematics.cxx:327
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
double Q2YtoX(double Ev, double M, double Q2, double y)
Definition: KineUtils.cxx:1117
const InitialState & InitState(void) const
Definition: Interaction.h:69
const ProcessInfo & ProcInfo(void) const
Definition: Interaction.h:70
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
const Target & Tgt(void) const
Definition: InitialState.h:67
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
bool IsCoherent(void) const
Definition: ProcessInfo.cxx:97
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.

1204 {
1205  Kinematics * kine = in->KinePtr();
1206 
1207  if(kine->KVSet(kKVW) && kine->KVSet(kKVQ2)) {
1208  const InitialState & init_state = in->InitState();
1209  double Ev = init_state.ProbeE(kRfHitNucRest);
1210  double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off mass shell
1211  double W = kine->W();
1212  double Q2 = kine->Q2();
1213 
1214  double x=-1,y=-1;
1215  kinematics::WQ2toXY(Ev,M,W,Q2,x,y);
1216  kine->Setx(x);
1217  kine->Sety(y);
1218  }
1219 }
double W(bool selected=false) const
Definition: Kinematics.cxx:167
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
Kinematics * KinePtr(void) const
Definition: Interaction.h:76
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
bool KVSet(KineVar_t kv) const
Definition: Kinematics.cxx:327
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
Definition: KineUtils.cxx:1046
void Setx(double x, bool selected=false)
Definition: Kinematics.cxx:241
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
void Sety(double y, bool selected=false)
Definition: Kinematics.cxx:253
const InitialState & InitState(void) const
Definition: Interaction.h:69
double Q2(bool selected=false) const
Definition: Kinematics.cxx:135
const Target & Tgt(void) const
Definition: InitialState.h:67
#define W(x)
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
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().

1016 {
1017  const ProcessInfo & process_info = interaction->ProcInfo();
1018 
1019  if(process_info.IsQuasiElastic()) {
1020  // hadronic inv. mass is equal to the recoil nucleon on-shell mass
1021  int rpdgc = interaction->RecoilNucleonPdg();
1022  double M = PDGLibrary::Instance()->Find(rpdgc)->Mass();
1023  return M;
1024  }
1025 
1026  const Kinematics & kinematics = interaction->Kine();
1027  if(kinematics.KVSet(kKVW)) {
1028  double W = kinematics.W();
1029  return W;
1030  }
1031  if(kinematics.KVSet(kKVx) && kinematics.KVSet(kKVy)) {
1032  const InitialState & init_state = interaction->InitState();
1033  double Ev = init_state.ProbeE(kRfHitNucRest);
1034  double M = init_state.Tgt().HitNucP4Ptr()->M();
1035  double M2 = M*M;
1036  double x = kinematics.x();
1037  double y = kinematics.y();
1038  double W2 = TMath::Max(0., M2 + 2*Ev*M*y*(1-x));
1039  double W = TMath::Sqrt(W2);
1040  return W;
1041  }
1042  SLOG("KineLimits", pERROR) << "Couldn't compute W for \n"<< *interaction;
1043  return 0;
1044 }
double W(bool selected=false) const
Definition: Kinematics.cxx:167
#define pERROR
Definition: Messenger.h:60
bool IsQuasiElastic(void) const
Definition: ProcessInfo.cxx:67
Generated/set kinematical variables for an event.
Definition: Kinematics.h:40
double x(bool selected=false) const
Definition: Kinematics.cxx:109
double y(bool selected=false) const
Definition: Kinematics.cxx:122
bool KVSet(KineVar_t kv) const
Definition: Kinematics.cxx:327
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
Definition: ProcessInfo.h:44
TLorentzVector * HitNucP4Ptr(void) const
Definition: Target.cxx:264
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
const Target & Tgt(void) const
Definition: InitialState.h:67
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Definition: Messenger.h:85
void kinematics()
Definition: kinematics.C:10
#define W(x)
double ProbeE(RefFrame_t rf) const
Initial State information.
Definition: InitialState.h:49
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().

1048 {
1049 // Converts (W,Q2) => (x,y)
1050 // Uses the system: a) W^2 - M^2 = 2*Ev*M*y*(1-x) and b) Q^2 = 2*x*y*M*Ev
1051 // Ev is the neutrino energy at the struck nucleon rest frame
1052 // M is the nucleon mass - it does not need to be on the mass shell
1053 
1054  double M2 = TMath::Power(M,2);
1055  double W2 = TMath::Power(W,2);
1056 
1057  x = Q2 / (W2-M2+Q2);
1058  y = (W2-M2+Q2) / (2*M*Ev);
1059 
1060  x = TMath::Min(1.,x);
1061  y = TMath::Min(1.,y);
1062  x = TMath::Max(0.,x);
1063  y = TMath::Max(0.,y);
1064 
1065  LOG("KineLimits", pDEBUG)
1066  << "(W=" << W << ",Q2=" << Q2 << ") => (x="<< x << ", y=" << y<< ")";
1067 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define W(x)
#define pDEBUG
Definition: Messenger.h:64
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().

1105 {
1106 // Converts (x,y) => Q2
1107 // Ev is the neutrino energy at the struck nucleon rest frame
1108 // M is the nucleon mass - it does not need to be on the mass shell
1109 
1110  double Q2 = 2*x*y*M*Ev;
1111 
1112  LOG("KineLimits", pDEBUG) << "(x=" << x << ",y=" << y << ") => Q2=" << Q2;
1113 
1114  return Q2;
1115 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
#define pDEBUG
Definition: Messenger.h:64
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().

1089 {
1090 // Converts (x,y) => W
1091 // Ev is the neutrino energy at the struck nucleon rest frame
1092 // M is the nucleon mass - it does not need to be on the mass shell
1093 
1094  double M2 = TMath::Power(M,2);
1095  double W2 = M2 + 2*Ev*M*y*(1-x);
1096  double W = TMath::Sqrt(TMath::Max(0., W2));
1097 
1098  LOG("KineLimits", pDEBUG) << "(x=" << x << ",y=" << y << ") => W=" << W;
1099 
1100  return W;
1101 }
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
#define W(x)
#define pDEBUG
Definition: Messenger.h:64
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().

1071 {
1072 // Converts (x,y) => (W,Q2)
1073 // Uses the system: a) W^2 - M^2 = 2*Ev*M*y*(1-x) and b) Q^2 = 2*x*y*M*Ev
1074 // Ev is the neutrino energy at the struck nucleon rest frame
1075 // M is the nucleon mass - it does not need to be on the mass shell
1076 
1077  double M2 = TMath::Power(M,2);
1078  double W2 = M2 + 2*Ev*M*y*(1-x);
1079 
1080  W = TMath::Sqrt(TMath::Max(0., W2));
1081  Q2 = 2*x*y*M*Ev;
1082 
1083  LOG("KineLimits", pDEBUG)
1084  << "(x=" << x << ",y=" << y << " => (W=" << W << ",Q2=" << Q2 << ")";
1085 }
double Q2(const Interaction *const i)
Definition: KineUtils.cxx:991
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Definition: Messenger.h:97
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Definition: Var.cxx:326
#define W(x)
#define pDEBUG
Definition: Messenger.h:64