30 using namespace genie;
69 double dW = (W.
max-W.
min)/(kNW-1);
71 for(
int iw=0; iw<kNW; iw++) {
74 double dQ2 = (Q2.
max-Q2.
min);
95 const double kdV = kdx*kdy;
97 double cW=-1, cQ2 = -1;
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;
121 <<
"Couldn't compute phase space volume for " 151 <<
"Computing Jacobian for transformation: " 182 J = 1. / (kine.
x() * kine.
y());
191 J = 1. / (kine.
Q2() * kine.
y());
234 J = TMath::Power(2*M*Ev,2) *
y;
248 J = 2*TMath::Power(M*Ev,2) * y/
W;
259 <<
"*** Can not compute Jacobian for transforming: " 266 if(!forward) J = 1./
J;
277 bool matched = ( (a==inpa&&b==inpb) || (a==inpb&&b==inpa) );
280 if(a==inpa&&b==inpb) fwd =
true;
292 double M2 = TMath::Power(M,2);
293 double s = M2 + 2*M*Ev;
310 double Ev,
double M,
double ml,
double W,
double Q2min_cut)
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;
327 double auxC = 0.5*(s-M2)/s;
328 double aux1 = s + ml2 - W2;
329 double aux2 = aux1*aux1 - 4*s*ml2;
331 (aux2 < 0) ? ( aux2 = 0 ) : ( aux2 =
TMath::Sqrt(aux2) );
333 Q2.
max = -ml2 + auxC * (aux1 + aux2);
334 Q2.
min = -ml2 + auxC * (aux1 - aux2);
337 Q2.
max = TMath::Max(0., Q2.
max);
338 Q2.
min = TMath::Max(0., Q2.
min);
341 if(Q2.
min < Q2min_cut) {Q2.
min = Q2min_cut; }
348 double Ev,
double M,
double ml,
double W,
double q2min_cut)
360 double Ev,
double M,
double ml,
double Q2min_cut)
369 if(W.
min<0)
return Q2;
376 double Ev,
double M,
double ml,
double q2min_cut)
392 double M2 = TMath::Power(M, 2.);
393 double ml2 = TMath::Power(ml,2.);
394 double s = M2 + 2*M*Ev;
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);
423 for(
unsigned int i=0;
i<N;
i++) {
424 double x = TMath::Power(10, logxmin +
i*dlogx);
431 if(y.
max >= 0 && y.
max <= 1 && y.
min >= 0 && y.
min <= 1) {
443 double Ev,
double M,
double ml,
double x)
452 double Ev2 = TMath::Power(Ev,2);
453 double ml2 = TMath::Power(ml,2);
461 double a = 0.5 * ml2/(M*Ev*
x);
463 double c = 1 + 0.5*x*M/Ev;
464 double d = TMath::Max(0., TMath::Power(1-a,2.) - b);
466 double A = 0.5 * (1-a-0.5*
b)/c;
481 double M2 = TMath::Power(M,2);
482 double ml2 = TMath::Power(ml,2);
483 double s = M2 + 2*M*El + ml2;
500 double El,
double ml,
double M,
double W)
508 double M2 = TMath::Power(M, 2.);
509 double ml2 = TMath::Power(ml, 2.);
510 double W2 = TMath::Power(W, 2.);
511 double s = M2 + 2*M*El + ml2;
517 double auxC = 0.5*(s - M2 - ml2)/s;
518 double aux1 = s + ml2 - W2;
519 double aux2 = aux1*aux1 - 4*s*ml2;
521 (aux2 < 0) ? ( aux2 = 0 ) : ( aux2 =
TMath::Sqrt(aux2) );
523 Q2.
max = -ml2 + auxC * (aux1 + aux2);
524 Q2.
min = -ml2 + auxC * (aux1 - aux2);
527 Q2.
max = TMath::Max(0., Q2.
max);
528 Q2.
min = TMath::Max(0., Q2.
min);
538 double El,
double ml,
double M,
double W)
550 double El,
double ml,
double M)
559 if(W.
min<0)
return Q2;
566 double El,
double ml,
double M)
582 double M2 = TMath::Power(M, 2.);
583 double ml2 = TMath::Power(ml,2.);
584 double s = M2 + 2*M*El + ml2;
608 const unsigned int N=100;
609 const double logxmin = TMath::Log10(xl.
min);
610 const double logxmax = TMath::Log10(xl.
max);
611 const double dlogx = (logxmax-logxmin) / (
double)(N-1);
613 for(
unsigned int i=0;
i<N;
i++) {
614 double x = TMath::Power(10, logxmin +
i*dlogx);
621 if(y.
max >= 0 && y.
max <= 1 && y.
min >= 0 && y.
min <= 1) {
633 double El,
double ml,
double M,
double x)
642 double El2 = TMath::Power(El,2);
643 double ml2 = TMath::Power(ml,2);
651 double a = 0.5 * ml2/(M*El*
x);
653 double c = 1 + 0.5*x*M/El;
654 double d = TMath::Max(0., TMath::Power(1-a,2.) - b);
656 double A = 0.5 * (1-a-0.5*
b)/c;
682 double Mn2 = Mn * Mn;
683 double mlep2 = mlep * mlep;
684 double s = Mn2 + 2.0 * Mn * Ev;
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;
693 double A = (s - Mn * Mn) / 2.0;
695 double C = 0.5 * (W2min + mlep2 - Mn2 * (W2min - mlep2) / s );
698 <<
"Q2 kinematic limits calculation failed for CohQ2Lim. " 699 <<
"Assuming Q2min = 0.0";
701 Q2.
min = TMath::Max(0., A * B - C);
704 <<
"Q2 kinematic limits calculation failed for CohQ2Lim. " 705 <<
"Assuming Q2min = 0.0";
721 double Ev,
double Q2)
730 double s = Mn * Mn + 2.0 * Mn * Ev;
731 double Mnterm = 1 - Mn * Mn /
s;
732 double Mlterm = 1 - mlep * mlep /
s;
734 double T1 = 0.25 * s * s * Mnterm * Mnterm * Mlterm;
735 double T2 = Q2 - (0.5 * s * Mnterm) + (0.5 * mlep * mlep * Mnterm);
738 W2l.
max = (T1 - T2 * T2 ) *
740 (1.0 / (Q2 + mlep * mlep));
746 double Q2,
double Mn,
double xsi)
752 double nu_min = (W2min + Q2 - Mn * Mn) / (2.0 * Mn);
753 double nu_max = (W2max + Q2 - Mn * Mn) / (2.0 * Mn);
756 nul.
min = (xsiQ > nu_min) ? xsiQ : nu_min;
763 double Ev,
double Q2,
double xsi)
770 if (W2lim.
min > W2lim.
max) {
772 <<
"Kinematically forbidden region in CohYLim. W2min = " << W2lim.
min 773 <<
"; W2max =" << W2lim.
max;
775 <<
" Mn = " << Mn <<
"; mpi = " << mpi <<
"; mlep = " 776 << mlep <<
"; Ev = " << Ev <<
"; Q2 = " <<
Q2;
781 ylim.
min = nulim.
min / Ev;
782 ylim.
max = nulim.
max / Ev;
803 return (Mn + mpi) * (Mn + mpi);
810 double M2 = TMath::Power(M,2);
811 double ml2 = TMath::Power(ml,2);
812 double s = M2 + 2*M*Ev + ml2;
830 double Ev,
double M,
double ml,
double W,
double Q2min_cut)
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;
844 double E1CM = (s + ml2 - M2) / (2.*sqs);
845 double p1CM = TMath::Max(0., E1CM*E1CM - ml2);
847 double E3CM = (s + ml2 - W2) / (2.*sqs);
848 double p3CM = TMath::Max(0., E3CM*E3CM - ml2);
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;
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);
863 SLOG(
"KineLimits",
pDEBUG) <<
"Nominal Q^2 limits: " << Q2.
min <<
" , " << Q2.
max;
865 Q2.
max = TMath::Max(0., Q2.
max);
866 Q2.
min = TMath::Max(0., Q2.
min);
869 if(Q2.
min < Q2min_cut) {Q2.
min = Q2min_cut; }
876 double Ev,
double M,
double ml,
double W,
double q2min_cut)
888 double Ev,
double M,
double ml,
double Q2min_cut)
897 if(W.
min<0)
return Q2;
904 double Ev,
double M,
double ml,
double q2min_cut)
922 double W2min = Wmin*
Wmin;
923 SLOG(
"KineLimits",
pDEBUG) <<
"W^2_min = " << W2min;
925 SLOG(
"KineLimits",
pDEBUG) <<
"Q^2 range : " << Q2l.
min <<
" , " << Q2l.
max;
928 x.
min = Q2l.
min / (Q2l.
min + W2min - M2);
929 x.
max = Q2l.
max / (Q2l.
max + W2min - M2);
940 double W2min = Wmin*
Wmin;
944 y.
min = (Q2l.
min + W2min - M2) / (2*Ev*M);
945 y.
max = (Q2l.
max + W2min - M2) / (2*Ev*M);
952 double Ev,
double M,
double ml,
double x)
964 double W2min = Wmin*
Wmin;
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);
999 double Q2 = kinematics.
Q2();
1005 double x = kinematics.
x();
1006 double y = kinematics.
y();
1008 double Q2 = 2*Mn*Ev*x*
y;
1011 SLOG(
"KineLimits",
pERROR) <<
"Couldn't compute Q^2 for \n"<< *interaction;
1028 double W = kinematics.
W();
1036 double x = kinematics.
x();
1037 double y = kinematics.
y();
1038 double W2 = TMath::Max(0., M2 + 2*Ev*M*y*(1-x));
1042 SLOG(
"KineLimits",
pERROR) <<
"Couldn't compute W for \n"<< *interaction;
1047 double Ev,
double M,
double W,
double Q2,
double &
x,
double &
y)
1054 double M2 = TMath::Power(M,2);
1055 double W2 = TMath::Power(W,2);
1057 x = Q2 / (W2-M2+
Q2);
1058 y = (W2-M2+
Q2) / (2*M*Ev);
1060 x = TMath::Min(1.,x);
1061 y = TMath::Min(1.,y);
1062 x = TMath::Max(0.,x);
1063 y = TMath::Max(0.,y);
1066 <<
"(W=" << W <<
",Q2=" << Q2 <<
") => (x="<< x <<
", y=" << y<<
")";
1070 double Ev,
double M,
double &
W,
double &
Q2,
double x,
double y)
1077 double M2 = TMath::Power(M,2);
1078 double W2 = M2 + 2*Ev*M*y*(1-
x);
1084 <<
"(x=" << x <<
",y=" << y <<
" => (W=" << W <<
",Q2=" << Q2 <<
")";
1088 double Ev,
double M,
double x,
double y)
1094 double M2 = TMath::Power(M,2);
1095 double W2 = M2 + 2*Ev*M*y*(1-
x);
1098 LOG(
"KineLimits",
pDEBUG) <<
"(x=" << x <<
",y=" << y <<
") => W=" <<
W;
1104 double Ev,
double M,
double x,
double y)
1110 double Q2 = 2*x*y*M*Ev;
1112 LOG(
"KineLimits",
pDEBUG) <<
"(x=" << x <<
",y=" << y <<
") => Q2=" <<
Q2;
1118 double Ev,
double M,
double Q2,
double y)
1123 assert(Ev > 0. && M > 0. && Q2 > 0. && y > 0.);
1125 double x = Q2 / (2. * y * M * Ev);
1127 LOG(
"KineLimits",
pDEBUG) <<
"(Ev=" << Ev <<
",Q2=" << Q2
1128 <<
",y=" << y <<
",M=" << M <<
") => x=" <<
x;
1134 double x,
double Q2,
double M,
double mc)
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);
1148 if(xc>=1 || W<=Wmin)
return false;
1153 double x,
double Q2,
double M,
double mc)
1160 double mc2 = TMath::Power(mc,2);
1161 double v = 0.5*Q2/(M*
x);
1162 double xc = x + 0.5*mc2/(M*
v);
1178 if (min_cut > range.
max || max_cut < range.
min) {
1193 double x = kine->
x();
1194 double y = kine->
y();
1211 double W = kine->
W();
1212 double Q2 = kine->
Q2();
1241 double y = kine->
y();
1242 double Q2 = kine->
Q2();
1250 double *
x,
double *
par)
1259 double xsmax = par[2];
1260 double Wmax = par[3];
1275 double hwfe = mD+2*gD;
1276 double lwfe = mD-gD/2;
1280 func = xsmax / (1 + 5* TMath::Power((W-lwfe)/gD,2));
1281 }
else if (W > hwfe) {
1283 func = xsmax / (1 + TMath::Power((W-hwfe)/gD,2));
1293 double plateau_edge = Wmax-0.1;
1295 if (W > plateau_edge) {
1300 func = xsmax / (1 + TMath::Power((W-plateau_edge)/gD,2));
1305 func *= (1 - (0.5-QD2));
1312 double *
x,
double *
par)
1319 double xpeak = par[0];
1320 double xmin = par[1];
1322 double xsmax = par[2];
1326 if(xb < xpeak/20.) {
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) {
1333 double plateau_edge = 2*xpeak;
1334 double slope = xsmax/(xmax - plateau_edge);
1335 func = xsmax - slope * (xb - plateau_edge);
1344 double *
x,
double *
par)
1351 double xsmax = 3*par[0];
1354 if(yb<0.|| yb>1.)
return 0.;
1355 if(xb<0.|| xb>1.)
return 0.;
1357 if(Ev<1)
return xsmax;
1358 if(xb/Ev<1E-4 && yb>0.95)
return 5*xsmax;
1362 double yp = (Ev>2.5) ? 2.5/Ev : 1;
1369 xs0 = xsmax - (yb-yp)*xsmax;
1371 double d = TMath::Power( (xb-xp)/0.075, 2);
1375 func = xsmax - (yb-yp)*xsmax;
double COHImportanceSamplingEnvelope(double *x, double *par)
static const double kMaxX
double W(bool selected=false) const
const KPhaseSpace & PhaseSpace(void) const
bool TransformMatched(KinePhaseSpace_t ia, KinePhaseSpace_t ib, KinePhaseSpace_t a, KinePhaseSpace_t b, bool &fwd)
std::map< std::string, double > xmax
double J(double q0, double q3, double Enu, double ml)
THE MAIN GENIE PROJECT NAMESPACE
Range1D_t InelXLim(double El, double ml, double M)
Range1D_t InelWLim(double El, double ml, double M)
void UpdateXYFromWQ2(const Interaction *in)
double Q2(const Interaction *const i)
Range1D_t Inelq2Lim_W(double El, double ml, double M, double W)
double DISImportanceSamplingEnvelope(double *x, double *par)
void SetQ2(double Q2, bool selected=false)
Kinematics * KinePtr(void) const
Range1D_t InelWLim(double Ev, double M, double ml)
int RecoilNucleonPdg(void) const
recoil nucleon pdg
A simple [min,max] interval for doubles.
Range1D_t Darkq2Lim(double Ev, double M, double ml, double q2min_cut=-1 *controls::kMinQ2Limit)
bool IsQuasiElastic(void) const
static const double kMaxY
Range1D_t DarkQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
bool IsWithinLimits(double x, Range1D_t range)
Generated/set kinematical variables for an event.
static const double kPhotontest
Range1D_t CohYLim(double Mn, double mpi, double mlep, double Ev, double Q2, double xsi)
double x(bool selected=false) const
static const double kLightestChmHad
static const double kMinQ2Limit
double CohW2Min(double Mn, double mpi)
static const double kMinY
Range1D_t Q2Lim_W(void) const
Q2 limits @ fixed W.
Range1D_t InelQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
enum genie::EKinePhaseSpace KinePhaseSpace_t
Range1D_t Limits(KineVar_t kvar) const
Return the kinematical variable limits.
double XYtoW(double Ev, double M, double x, double y)
double PhaseSpaceVolume(const Interaction *const i, KinePhaseSpace_t ps)
Range1D_t InelQ2Lim_W(double Ev, double M, double ml, double W, double Q2min_cut=controls::kMinQ2Limit)
double y(bool selected=false) const
double W(const Interaction *const i)
Summary information for an interaction.
double SlowRescalingVar(double x, double Q2, double M, double mc)
Range1D_t InelXLim(double Ev, double M, double ml)
Range1D_t DarkQ2Lim(double Ev, double M, double ml, double Q2min_cut=controls::kMinQ2Limit)
Range1D_t Inelq2Lim_W(double Ev, double M, double ml, double W, double q2min_cut=-1 *controls::kMinQ2Limit)
#define LOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a string (using the FILE...
Range1D_t CohNuLim(double W2min, double W2max, double Q2, double Mn, double xsi)
bool KVSet(KineVar_t kv) const
double QD2toQ2(double QD2)
Range1D_t Darkq2Lim_W(double Ev, double M, double ml, double W, double q2min_cut=-1 *controls::kMinQ2Limit)
static string AsString(KinePhaseSpace_t kps)
A class encapsulating an enumeration of interaction types (EM, Weak-CC, Weak-NC) and scattering types...
double XYtoQ2(double Ev, double M, double x, double y)
void XYtoWQ2(double Ev, double M, double &W, double &Q2, double x, double y)
void WQ2toXY(double Ev, double M, double W, double Q2, double &x, double &y)
const Kinematics & Kine(void) const
static const double kNeutronMass
Range1D_t Cohq2Lim(double Mn, double mpi, double mlep, double Ev)
Range1D_t CohQ2Lim(double Mn, double mpi, double mlep, double Ev)
static const double kMinX
double func(double x, double y)
static const double kASmallNum
bool IsAboveCharmThreshold(double x, double Q2, double M, double mc)
double Q2toQD2(double Q2)
Range1D_t InelQ2Lim(double El, double ml, double M)
Range1D_t Inelq2Lim(double El, double ml, double M)
void Setx(double x, bool selected=false)
static const double kMQD2
void UpdateWQ2FromXY(const Interaction *in)
Range1D_t InelYLim(double El, double ml, double M)
void SetW(double W, bool selected=false)
TLorentzVector * HitNucP4Ptr(void) const
Range1D_t InelQ2Lim_W(double El, double ml, double M, double W)
static PDGLibrary * Instance(void)
static const double kPionMass
Range1D_t InelYLim_X(double Ev, double M, double ml, double x)
double Q2YtoX(double Ev, double M, double Q2, double y)
Var Sqrt(const Var &v)
Use to take sqrt of a var.
Range1D_t DarkYLim(double Ev, double M, double ml)
void Sety(double y, bool selected=false)
void UpdateXFromQ2Y(const Interaction *in)
TParticlePDG * Find(int pdgc)
Range1D_t InelYLim_X(double El, double ml, double M, double x)
assert(nhit_max >=nhit_nbins)
double Jacobian(const Interaction *const i, KinePhaseSpace_t f, KinePhaseSpace_t t)
const InitialState & InitState(void) const
const ProcessInfo & ProcInfo(void) const
Range1D_t DarkYLim_X(double Ev, double M, double ml, double x)
Range1D_t Inelq2Lim(double Ev, double M, double ml, double q2min_cut=-1 *controls::kMinQ2Limit)
double Q2(bool selected=false) const
const Target & Tgt(void) const
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
#define SLOG(stream, priority)
A macro that returns the requested log4cpp::Category appending a short string (using the FUNCTION and...
Range1D_t InelYLim(double Ev, double M, double ml)
double ProbeE(RefFrame_t rf) const
Range1D_t DarkXLim(double Ev, double M, double ml)
Most commonly used PDG codes. A set of utility functions to handle PDG codes is provided in PDGUtils...
Range1D_t CohW2Lim(double Mn, double mpi, double mlep, double Ev, double Q2)
void ApplyCutsToKineLimits(Range1D_t &r, double min, double max)
Range1D_t DarkWLim(double Ev, double M, double ml)
Initial State information.
bool IsCoherent(void) const
double RESImportanceSamplingEnvelope(double *x, double *par)