Skip to content

Commit

Permalink
Fixes NuWro#11
Browse files Browse the repository at this point in the history
  • Loading branch information
@karuboniru
karuboniru committed Nov 28, 2023
1 parent 4d2d5a8 commit 9de466d
Showing 1 changed file with 37 additions and 32 deletions.
69 changes: 37 additions & 32 deletions src/dis/dis_nc.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,11 @@
//Cross section for inelastic scattering on proton
// nu + p --> nu +X (NC, lepton)

double xx_boundary_normalize(double xx_in){
if (xx_in <= 0) return nextafter(0., 1.);
if (xx_in >= 1) return nextafter(1., 0.);
return xx_in;
}

double
sigma_xy_nc_nu_p_Paschos (double E, double xx, double yy, double m)
Expand All @@ -36,7 +41,7 @@ sigma_xy_nc_nu_p_Paschos (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_Paschos (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -66,7 +71,7 @@ sigma_xy_nc_nu_p_grv (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p_grv (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_grv (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -96,7 +101,7 @@ sigma_xy_nc_nu_p_u (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p_u (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_u (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -127,7 +132,7 @@ sigma_xy_nc_nu_p_d (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p_d (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_d (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -161,7 +166,7 @@ sigma_xy_nc_nu_p_ubar (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p_ubar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_ubar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -196,7 +201,7 @@ sigma_xy_nc_nu_p_dbar (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p_dbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_dbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -227,7 +232,7 @@ sigma_xy_nc_nu_p_s (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p_s (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_s (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -262,7 +267,7 @@ sigma_xy_nc_nu_p_sbar (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_p_sbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_p_sbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -294,7 +299,7 @@ sigma_xy_nc_anu_p_Paschos (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_Paschos (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -326,7 +331,7 @@ sigma_xy_nc_anu_p_grv (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p_grv (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_grv (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand All @@ -353,7 +358,7 @@ sigma_xy_nc_anu_p_u (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p_u (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_u (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -381,7 +386,7 @@ sigma_xy_nc_anu_p_d (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p_d (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_d (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -409,7 +414,7 @@ sigma_xy_nc_anu_p_ubar (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p_ubar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_ubar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -438,7 +443,7 @@ sigma_xy_nc_anu_p_dbar (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p_dbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_dbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -466,7 +471,7 @@ sigma_xy_nc_anu_p_s (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p_s (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_s (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -495,7 +500,7 @@ sigma_xy_nc_anu_p_sbar (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_p_sbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_p_sbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -533,7 +538,7 @@ sigma_xy_nc_nu_n_Paschos (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_n (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_Paschos (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -584,7 +589,7 @@ sigma_xy_nc_nu_n_grv (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_n_grv (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_grv (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand All @@ -605,7 +610,7 @@ cr_sec_nc_nu_n_grv (double E, double W, double nu, double m)
double
cr_sec_nc_nu_n_u (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_u (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -636,7 +641,7 @@ sigma_xy_nc_nu_n_d (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_n_d (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_d (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -670,7 +675,7 @@ sigma_xy_nc_nu_n_ubar (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_n_ubar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_ubar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -705,7 +710,7 @@ sigma_xy_nc_nu_n_dbar (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_n_dbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_dbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -736,7 +741,7 @@ sigma_xy_nc_nu_n_s (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_n_s (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_s (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -771,7 +776,7 @@ sigma_xy_nc_nu_n_sbar (double E, double xx, double yy, double m)
double
cr_sec_nc_nu_n_sbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_nu_n_sbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -803,7 +808,7 @@ sigma_xy_nc_anu_n_Paschos (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_Paschos (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -835,7 +840,7 @@ sigma_xy_nc_anu_n_grv (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n_grv (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_grv (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand All @@ -862,7 +867,7 @@ sigma_xy_nc_anu_n_u (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n_u (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_u (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -890,7 +895,7 @@ sigma_xy_nc_anu_n_d (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n_d (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_d (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -918,7 +923,7 @@ sigma_xy_nc_anu_n_ubar (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n_ubar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_ubar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -947,7 +952,7 @@ sigma_xy_nc_anu_n_dbar (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n_dbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_dbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -975,7 +980,7 @@ sigma_xy_nc_anu_n_s (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n_s (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_s (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down Expand Up @@ -1004,7 +1009,7 @@ sigma_xy_nc_anu_n_sbar (double E, double xx, double yy, double m)
double
cr_sec_nc_anu_n_sbar (double E, double W, double nu, double m)
{
double xx = (M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu);
double xx = xx_boundary_normalize((M2 - W * W + 2 * M12 * nu) / (2 * M12 * nu));
double yy = nu / E;
return sigma_xy_nc_anu_n_sbar (E, xx, yy, m) * W / (M12 * E * nu) //jakobian
;
Expand Down

0 comments on commit 9de466d

Please sign in to comment.