diff --git a/reactions/AFB/include/ReactionAFB.h b/reactions/AFB/include/ReactionAFB.h index 7b9a74b0f72d8bc1f6b813ea6fd3f5955b859b92..b28e10f6a60983838fb92e8c1a18d59fee38ee8f 100644 --- a/reactions/AFB/include/ReactionAFB.h +++ b/reactions/AFB/include/ReactionAFB.h @@ -38,7 +38,12 @@ class ReactionAFB : public ReactionTheory ReactionTheory* ptr; }; + static string integration_param; + static int integration_switch; + static int key_param; + static size_t alloc_space; static size_t calls; + static double epsabs; static double epsrel; static double PI; diff --git a/reactions/AFB/src/ReactionAFB.cc b/reactions/AFB/src/ReactionAFB.cc index a8808fe489dad28acbdd46d44b97de174e4f888e..6b3c93c9948222dd30a35a7c9edd6f582617f346 100644 --- a/reactions/AFB/src/ReactionAFB.cc +++ b/reactions/AFB/src/ReactionAFB.cc @@ -8,6 +8,7 @@ #include "ReactionAFB.h" #include "iostream" +#include "cstring" #include <gsl/gsl_integration.h> using namespace std; @@ -17,6 +18,10 @@ double ReactionAFB::PI; double ReactionAFB::GeVtofb_param, ReactionAFB::alphaEM_param, ReactionAFB::stheta2W_param, ReactionAFB::MZ_param, ReactionAFB::GammaZ_param; double ReactionAFB::energy_param, ReactionAFB::eta_cut_param, ReactionAFB::pT_cut_param, ReactionAFB::y_min_param, ReactionAFB::y_max_param; +int ReactionAFB::integration_switch; +string ReactionAFB::integration_param; +int ReactionAFB::key_param; + double ReactionAFB::e_param, ReactionAFB::gsm_param, ReactionAFB::smangle_param; double ReactionAFB::foton_Vu, ReactionAFB::foton_Au, ReactionAFB::foton_Vd, ReactionAFB::foton_Ad, ReactionAFB::foton_Vl, ReactionAFB::foton_Al, ReactionAFB::foton_Vnu, ReactionAFB::foton_Anu; double ReactionAFB::Z_Vu, ReactionAFB::Z_Au, ReactionAFB::Z_Vd, ReactionAFB::Z_Ad, ReactionAFB::Z_Vl, ReactionAFB::Z_Al, ReactionAFB::Z_Vnu, ReactionAFB::Z_Anu; @@ -27,7 +32,7 @@ double ReactionAFB::epsabs = 0; double ReactionAFB::epsrel = 1e-2; size_t ReactionAFB::calls; - +size_t ReactionAFB::alloc_space = 1000; //// Function returning the combination of propagators double *ReactionAFB::propagators (double Minv) @@ -102,21 +107,31 @@ double ReactionAFB::integration_uubarEF_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::uubarEF_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::uubarEF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return 2*result; } @@ -125,15 +140,24 @@ double ReactionAFB::integration_uubarEF_y (double Minv, void * ptr) { double ReactionAFB::integration_uubarEF (double Minv_inf, double Minv_sup, void* ptr) { double result, error; - - gsl_function F; - F.function = &(ReactionAFB::integration_uubarEF_y); - F.params = ptr; - double inf = Minv_inf; double sup = Minv_sup; - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::integration_uubarEF_y); + F.params = ptr; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return result; } @@ -193,39 +217,58 @@ double ReactionAFB::integration_uubarEB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::uubarEB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - + gsl_function F; + F.function = &(ReactionAFB::uubarEB_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////UUBAR EVEN BACKWARD Integration in invariant mass double ReactionAFB::integration_uubarEB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; - + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; + gsl_function F; F.function = &(ReactionAFB::integration_uubarEB_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -285,22 +328,32 @@ double ReactionAFB::integration_uubarOF_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::uubarOF_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - + gsl_function F; + F.function = &(ReactionAFB::uubarOF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } @@ -308,15 +361,24 @@ double ReactionAFB::integration_uubarOF_y (double Minv, void * ptr) { double ReactionAFB::integration_uubarOF (double Minv_inf, double Minv_sup, void* ptr) { double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_uubarOF_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return result; } @@ -377,38 +439,57 @@ double ReactionAFB::integration_uubarOB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::uubarOB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::uubarOB_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////UUBAR ODD BACKWARD Integration in invariant mass double ReactionAFB::integration_uubarOB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_uubarOB_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return result; } @@ -470,20 +551,30 @@ double ReactionAFB::integration_ubaruEF_y (double Minv, void * ptr) { integrationParams.ptr = (ReactionTheory*) ptr; double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ubaruEF_funct); - F.params = &integrationParams; - double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::ubaruEF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return 2*result; } @@ -491,17 +582,26 @@ double ReactionAFB::integration_ubaruEF_y (double Minv, void * ptr) { ////UBARU EVEN FORWARD Integration in invariant mass double ReactionAFB::integration_ubaruEF (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_ubaruEF_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -561,21 +661,31 @@ double ReactionAFB::integration_ubaruEB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ubaruEB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } + + gsl_function F; + F.function = &(ReactionAFB::ubaruEB_funct); + F.params = &integrationParams; - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return 2*result; } @@ -583,16 +693,25 @@ double ReactionAFB::integration_ubaruEB_y (double Minv, void * ptr) { ////UBARU EVEN BACKWARD Integration in invariant mass double ReactionAFB::integration_ubaruEB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_ubaruEB_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return result; } @@ -653,21 +772,31 @@ double ReactionAFB::integration_ubaruOF_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ubaruOF_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::ubaruOF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return 2*result; } @@ -675,17 +804,26 @@ double ReactionAFB::integration_ubaruOF_y (double Minv, void * ptr) { ////UBARU ODD FORWARD Integration in invariant mass double ReactionAFB::integration_ubaruOF (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_ubaruOF_y); - F.params = ptr; + F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -745,39 +883,58 @@ double ReactionAFB::integration_ubaruOB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ubaruOB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::ubaruOB_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////UBARU ODD BACKWARD Integration in invariant mass double ReactionAFB::integration_ubaruOB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_ubaruOB_y); - F.params = ptr; + F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -839,39 +996,58 @@ double ReactionAFB::integration_ddbarEF_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ddbarEF_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::ddbarEF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////DDBAR EVEN FORWARD Integration in invariant mass double ReactionAFB::integration_ddbarEF (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; - + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; + gsl_function F; F.function = &(ReactionAFB::integration_ddbarEF_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -933,21 +1109,31 @@ double ReactionAFB::integration_ddbarEB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ddbarEB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } + + gsl_function F; + F.function = &(ReactionAFB::ddbarEB_funct); + F.params = &integrationParams; - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return 2*result; } @@ -955,17 +1141,26 @@ double ReactionAFB::integration_ddbarEB_y (double Minv, void * ptr) { ////DDBAR EVEN BACKWARD Integration in invariant mass double ReactionAFB::integration_ddbarEB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_ddbarEB_y); - F.params = ptr; + F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -1027,39 +1222,58 @@ double ReactionAFB::integration_ddbarOF_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ddbarOF_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::ddbarOF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////DDBAR ODD FORWARD Integration in invariant mass double ReactionAFB::integration_ddbarOF (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_ddbarOF_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -1121,39 +1335,58 @@ double ReactionAFB::integration_ddbarOB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::ddbarOB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::ddbarOB_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////DDBAR ODD BACKWARD Integration in invariant mass double ReactionAFB::integration_ddbarOB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_ddbarOB_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -1215,39 +1448,58 @@ double ReactionAFB::integration_dbardEF_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::dbardEF_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::dbardEF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////DBARD EVEN FORWARD Integration in invariant mass double ReactionAFB::integration_dbardEF (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_dbardEF_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return result; } @@ -1309,39 +1561,58 @@ double ReactionAFB::integration_dbardEB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::dbardEB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::dbardEB_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } ////DBARD EVEN BACKWARD Integration in invariant mass double ReactionAFB::integration_dbardEB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_dbardEB_y); F.params = ptr; - - double inf = Minv_inf; - double sup = Minv_sup; - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return result; } @@ -1403,22 +1674,32 @@ double ReactionAFB::integration_dbardOF_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::dbardOF_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::dbardOF_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } + return 2*result; } @@ -1426,16 +1707,25 @@ double ReactionAFB::integration_dbardOF_y (double Minv, void * ptr) { double ReactionAFB::integration_dbardOF (double Minv_inf, double Minv_sup, void* ptr) { double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_dbardOF_y); - F.params = ptr; + F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -1497,21 +1787,31 @@ double ReactionAFB::integration_dbardOB_y (double Minv, void * ptr) { integrationParams.Minv = Minv; integrationParams.ptr = (ReactionTheory*) ptr; - double result, error; - - gsl_function F; - F.function = &(ReactionAFB::dbardOB_funct); - F.params = &integrationParams; - + double result, error; double inf = y_min_param / log(energy_param/Minv); double sup; + if (y_max_param == 0.0) { sup = 1; } else { sup = y_max_param / log(energy_param/Minv); } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + gsl_function F; + F.function = &(ReactionAFB::dbardOB_funct); + F.params = &integrationParams; + + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } return 2*result; } @@ -1519,17 +1819,26 @@ double ReactionAFB::integration_dbardOB_y (double Minv, void * ptr) { ////DBARD ODD BACKWARD Integration in invariant mass double ReactionAFB::integration_dbardOB (double Minv_inf, double Minv_sup, void* ptr) { - double result, error; + double result, error; + double inf = Minv_inf; + double sup = Minv_sup; gsl_function F; F.function = &(ReactionAFB::integration_dbardOB_y); F.params = ptr; - double inf = Minv_inf; - double sup = Minv_sup; + if (integration_switch == 1) { + gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); + } + else if (integration_switch == 2) { + gsl_integration_workspace * w = gsl_integration_workspace_alloc(alloc_space); + gsl_integration_qag (&F, inf, sup, epsabs, epsrel, alloc_space, key_param, w, &result, &error); + gsl_integration_workspace_free (w); + } + else { + result = 0.0; + } - gsl_integration_qng (&F, inf, sup, epsabs, epsrel, &result, &error, &calls); - return result; } @@ -1574,6 +1883,7 @@ extern "C" ReactionAFB* create() { int ReactionAFB::initAtStart(const string &s) { // Parameters from "/reactions/AFB/yaml/parameters.yaml" + // Check energy parameter: std::cout << checkParam("energy") << std::endl; if ( ! checkParam("energy") ) { @@ -1606,6 +1916,32 @@ int ReactionAFB::initAtStart(const string &s) return 1; } + // Check integration routine parameters: + std::cout << checkParam("integration") << std::endl; + if ( ! checkParam("integration") ) { + std::cout << "\n\n FATAL ERROR: integration routine (integration) is not defined !!! \n\n" <<std::endl; + return 1; + } + integration_param = GetParamS("integration"); + if (not(integration_param.compare("QNG"))) { + integration_switch = 1; + } + else if (not(integration_param.compare("QAG"))) { + integration_switch = 2; + if ( ! checkParam("key") ) { + std::cout << "\n\n FATAL ERROR: rule for QAG integration is not defined !!! \n\n" <<std::endl; + return 1; + } + key_param = GetParamI("key"); + if ((key_param < 1) or (key_param > 6)) { + std::cout << "\n\n FATAL ERROR: rule for QAG integration has to be between 1 and 6 (1 <= key <= 6) !!! \n\n" <<std::endl; + return 1; + } + } + else { + std::cout << "\n\n FATAL ERROR: integration routine supported are QNG and QAG. Please select one of these two options !!! \n\n" <<std::endl; + } + // Constant PI = 3.14159265; @@ -1620,8 +1956,8 @@ int ReactionAFB::initAtStart(const string &s) energy_param = GetParam("energy"); eta_cut_param = GetParam("eta_cut"); pT_cut_param = GetParam("pT_cut"); - y_min_param = GetParam("y_min"); // not implemented yet - y_max_param = GetParam("y_max"); // not implemented yet + y_min_param = GetParam("y_min"); + y_max_param = GetParam("y_max"); // Calculate fixed parameters e_param = sqrt(4*PI*alphaEM_param); @@ -1681,12 +2017,18 @@ int ReactionAFB::compute(int dataSetID, valarray<double> &val, map<string, valar int Npnt_min = min.size(); int Npnt_max = max.size(); + + // check on the rapidity cut + if (y_min_param / log(energy_param/max[Npnt_max-1]) > 1) { + std::cout << "\n\nThe chosen lower rapidity cut is too high in this invariant mass range." << std::endl; + return 1; + } if (Npnt_min != Npnt_max) { std::cout << "\n\nFATAL ERROR: uneven number of Invariant mass min and max !!!" << std::endl; std::cout << "CHECK THE DATAFILE !!!" << std::endl; return 1; - } + } // Fill the array "val[i]" with the result of the AFB function for (int i = 0; i < Npnt_min; i++) { diff --git a/reactions/AFB/yaml/parameters.yaml b/reactions/AFB/yaml/parameters.yaml index c33fa29111d2cbf60f1f36ff36e8ed5b509ce0f0..55f43d90276175085f469b2743e1da0a93f67a86 100644 --- a/reactions/AFB/yaml/parameters.yaml +++ b/reactions/AFB/yaml/parameters.yaml @@ -1,7 +1,15 @@ ### Parameters for DY AFB +## Kinematic parameters energy : 13000.0 eta_cut : 2.5 pT_cut : 20.0 y_min : 0.0 y_max : 0.0 # if set to 0 no upper y_cut + + +## Integration routine +integration : "QNG" +# "QNG" for non-adaptive Gauss-Kronrod integration (faster, less precise); "QAG" for adaptive Gauss-Kronrod integration (slower, more precise) +key : 6 +# sets the precision for QAG integration. "1, 2, 3, 4, 5, 6" for increasing precision (corresponds to the 15, 21, 31, 41, 51, 61 point Gauss-Kronrod rule) diff --git a/reactions/APPLgrid/include/appl_grid/appl_grid.h b/reactions/APPLgrid/include/appl_grid/appl_grid.h deleted file mode 100644 index bc74a2050b09d2d0caa53b8fb56e4dd69b879f9e..0000000000000000000000000000000000000000 --- a/reactions/APPLgrid/include/appl_grid/appl_grid.h +++ /dev/null @@ -1,696 +0,0 @@ -// emacs: this is -*- c++ -*- - -// appl_grid.h - -// grid class header - all the functions needed to create and -// fill the grid from an NLO calculation program -// -// Copyright (C) 2007 Mark Sutton (sutt@hep.ucl.ac.uk) - -// $Id: appl_grid.h, v1.00 2007/10/16 17:01:39 sutt - -// Fixme: this needs to be tidied up. eg there are too many different, -// and too many version of, accessors for x/y, Q2/tau etc there -// should be only one set, for x and Q2 *or* y and tau, but -// not both. In fact they should be for x and Q2, since y and tau -// should be purely an internal grid issue of no concern for the -// user. - -#ifndef __APPL_GRID_H -#define __APPL_GRID_H - -#include <vector> -#include <iostream> -#include <sstream> -#include <cmath> -#include <string> -#include <exception> - -#include "TH1D.h" - - -double _fy(double x); -double _fx(double y); -double _fun(double y); - - -#include "correction.h" - -namespace appl { - - -/// forward declarations - full definitions included -/// from appl_grid.cxx -class igrid; -class appl_pdf; - - -const int MAXGRIDS = 5; - - -/// externally visible grid class -class grid { - -public: - - // grid error exception - class exception : public std::exception { - public: - exception(const std::string& s) { std::cerr << what() << " " << s << std::endl; }; - //exception(std::ostream& s) { std::cerr << what() << " " << s << std::endl; }; - exception(std::ostream& s) { std::stringstream ss; ss << s.rdbuf(); std::cerr << what() << " " << ss.str() << std::endl; }; - virtual const char* what() const throw() { return "appl::grid::exception"; } - }; - - typedef enum { STANDARD=0, AMCATNLO=1, SHERPA=2, LAST_TYPE=3 } CALCULATION; - -public: - - grid(int NQ2=50, double Q2min=10000.0, double Q2max=25000000.0, int Q2order=5, - int Nx=50, double xmin=1e-5, double xmax=0.9, int xorder=5, - int Nobs=20, double obsmin=100.0, double obsmax=7000.0, - std::string genpdf="mcfm_pdf", - int leading_order=0, int nloops=1, - std::string transform="f2"); - - grid( int Nobs, const double* obsbins, - int NQ2=50, double Q2min=10000.0, double Q2max=25000000.0, int Q2order=5, - int Nx=50, double xmin=1e-5, double xmax=0.9, int xorder=5, - std::string genpdf="mcfm_pdf", - int leading_order=0, int nloops=1, - std::string transform="f2" ); - - grid( const std::vector<double>& obs, - int NQ2=50, double Q2min=10000.0, double Q2max=25000000.0, int Q2order=5, - int Nx=50, double xmin=1e-5, double xmax=0.9, int xorder=5, - std::string genpdf="mcfm_pdf", - int leading_order=0, int nloops=1, - std::string transform="f2" ); - - // build a grid but don't build the internal igrids - these can be added later - grid( const std::vector<double>& obs, - std::string genpdf="nlojet_pdf", - int leading_order=0, int nloops=1, - std::string transform="f2" ); - - // copy constructor - grid(const grid& g); - - // read from a file - grid(const std::string& filename="./grid.root", const std::string& dirname="grid"); - - // add an igrid for a given bin and a given order - void add_igrid(int bin, int order, igrid* g); - - virtual ~grid(); - - // update grid with one set of event weights - void fill(const double x1, const double x2, const double Q2, - const double obs, - const double* weight, const int iorder); - - - void fill_phasespace(const double x1, const double x2, const double Q2, - const double obs, - const double* weight, const int iorder); - - - void fill_grid(const double x1, const double x2, const double Q2, - const double obs, - const double* weight, const int iorder) { - if (isOptimised()) fill(x1, x2, Q2, obs, weight, iorder); - else fill_phasespace(x1, x2, Q2, obs, weight, iorder); - } - - - void fill_index(const int ix1, const int ix2, const int iQ2, - const int iobs, - const double* weight, const int iorder); - - - // trim/untrim the grid to reduce memory footprint - void trim(); - void untrim(); - - // formatted output - std::ostream& print(std::ostream& s=std::cout) const; - - // don't do anything anymore - // void setuppdf(void (*pdf)(const double& , const double&, double* ) ); - - // get the interpolated pdf's - // void pdfinterp(double x1, double Q2, double* f); - - - // perform the convolution to a specified number of loops - // nloops=-1 gives the nlo part only - std::vector<double> vconvolute(void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, - double fscale_factor=1, - double Escale=1 ); - - std::vector<double> vconvolute(void (*pdf1)(const double& , const double&, double* ), - void (*pdf2)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, - double fscale_factor=1, - double Escale=1 ); - - - // perform the convolution to a specified number of loops - // nloops=-1 gives the nlo part only - std::vector<double> vconvolute(double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, - double fscale_factor=1 ) { - return vconvolute(pdf, alphas, nloops, rscale_factor, fscale_factor, Escale); - } - - - // perform the convolution to the max number of loops in grid - std::vector<double> vconvolute(void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return vconvolute( pdf, alphas, m_order-1 ); - } - - - // perform the convolution to the max number of loops in grid - std::vector<double> vconvolute(double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return vconvolute( Escale, pdf, alphas, m_order-1 ); - } - - - // perform the convolution to a specified number of loops - // for a single sub process, nloops=-1 gives the nlo part only - std::vector<double> vconvolute_subproc(int subproc, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, double Escale=1 ); - - - // perform the convolution to a specified number of loops - // for a single sub process, nloops=-1 gives the nlo part only - std::vector<double> vconvolute_subproc(int subproc, - double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1 ) { - return vconvolute_subproc(subproc, pdf, alphas, nloops, rscale_factor, Escale); - } - - - // perform the convolution to the max number of loops in grid - // for a single sub process - std::vector<double> vconvolute_subproc(int subproc, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return vconvolute_subproc( subproc, pdf, alphas, m_order-1 ); - } - - // perform the convolution to the max number of loops in grid - // for a single sub process - std::vector<double> vconvolute_subproc(int subproc, - double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return vconvolute_subproc( subproc, Escale, pdf, alphas, m_order-1 ); - } - - - double vconvolute_bin( int bin, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double&) ); - - - // perform the convolution to a specified number of loops - // nloops=-1 gives the nlo part only - TH1D* convolute(void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, - double fscale_factor=1, - double Escale=1 ); - - // perform the convolution to a specified number of loops - // nloops=-1 gives the nlo part only - TH1D* convolute(void (*pdf1)(const double& , const double&, double* ), - void (*pdf2)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, - double fscale_factor=1, - double Escale=1 ); - - - TH1D* convolute(double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, - double fscale_factor=1 ) { - return convolute(pdf, alphas, nloops, rscale_factor, fscale_factor, Escale); - } - - - // perform the convolution to the max number of loops in grid - TH1D* convolute(void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return convolute( pdf, alphas, m_order-1 ); - } - - // perform the convolution to the max number of loops in grid - TH1D* convolute(double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return convolute( Escale, pdf, alphas, m_order-1 ); - } - - - // perform the convolution to a specified number of loops - // for a single sub process, nloops=-1 gives the nlo part only - TH1D* convolute_subproc(int subproc, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1, double Escale=1 ); - - TH1D* convolute_subproc(int subproc, - double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ), - int nloops, - double rscale_factor=1 ) { - return convolute_subproc( subproc, pdf, alphas, nloops, rscale_factor, Escale); - } - - // perform the convolution to the max number of loops in grid - // for a single sub process - TH1D* convolute_subproc(int subproc, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return convolute_subproc( subproc, pdf, alphas, m_order-1 ); - } - - TH1D* convolute_subproc(int subproc, - double Escale, - void (*pdf)(const double& , const double&, double* ), - double (*alphas)(const double& ) ) { - return convolute_subproc( subproc, Escale, pdf, alphas, m_order-1 ); - } - - - // optimise the bin limits - void optimise(bool force=false); - void optimise(int NQ2, int Nx); - void optimise(int NQ2, int Nx1, int Nx2); - - // redefine the limits by hand - void redefine(int iobs, int iorder, - int NQ2, double Q2min, double Q2max, - int Nx, double xmin, double xmax); - - bool setNormalised(bool t=true) { return m_normalised=t; } - bool getNormalised() const { return m_normalised; } - - - // set the filling to be symmetric and test status - bool symmetrise(bool t=true) { return m_symmetrise=t; } - bool isSymmetric() const { return m_symmetrise; } - - bool reweight(bool t=false); - - // access to internal grids if need be - const igrid* weightgrid(int iorder, int iobs) const { return m_grids[iorder][iobs]; } - - // save grid to specified file - void Write(const std::string& filename, const std::string& dirname="grid", const std::string& pdfname="" ); - - // accessors for the observable after possible bin combination - int Nobs() const { return m_obs_bins_combined->GetNbinsX(); } - double obs(int iobs) const { return m_obs_bins_combined->GetBinCenter(iobs+1); } - int obsbin(double obs) const { return m_obs_bins_combined->FindBin(obs)-1; } - double obslow(int iobs) const { return m_obs_bins_combined->GetBinLowEdge(iobs+1); } - double obsmin() const { return obslow(0); } - double obsmax() const { return obslow(Nobs()); } - double deltaobs(int iobs) const { return m_obs_bins_combined->GetBinWidth(iobs+1); } - - const TH1D* getReference() const { return m_obs_bins_combined; } - TH1D* getReference() { return m_obs_bins_combined; } - - - // TH1D* getXReference() { - // combineReference(); - // return m_obs_bins_combined; - // } - - - // accessors for the observable befor any bin combination - int Nobs_internal() const { return m_obs_bins->GetNbinsX(); } - double obs_internal(int iobs) const { return m_obs_bins->GetBinCenter(iobs+1); } - int obsbin_internal(double obs) const { return m_obs_bins->FindBin(obs)-1; } - double obslow_internal(int iobs) const { return m_obs_bins->GetBinLowEdge(iobs+1); } - double deltaobs_internal(int iobs) const { return m_obs_bins->GetBinWidth(iobs+1); } - double obsmin_internal() const { return obslow_internal(0); } - double obsmax_internal() const { return obslow_internal(Nobs_internal()); } - - const TH1D* getReference_internal() const { return m_obs_bins; } - TH1D* getReference_internal() { return m_obs_bins; } - - - - - // number of subprocesses - int subProcesses(int i) const; - - // general status accessors - double& run() { return m_run; } - - // accessors for the status information - bool isOptimised() const { return m_optimised; } - bool isTrimmed() const { return m_trimmed; } - - // lowest order of process - int leadingOrder() const { return m_leading_order; } - - /// maximum number of orders ( lo=1, nlo=2, nnlo=3 ) - /// but aMC@NLO uses 4 grids for the NLO, so m_order - /// will be 4, but really it is still only available - /// 1 loop, so take account of this - int nloops() const { - if ( m_type!=AMCATNLO ) return m_order-1; - else if ( m_order>0 ) return 1; - else return 0; - } - - // find out which transform and which pdf combination are being used - std::string getTransform() const { return m_transform; } - - static double transformvar(); - static double transformvar(double v); - - std::string getGenpdf() const { return m_genpdfname; } - - std::string version() const { return m_version; } - std::string appl_version() const; - - double getCMSScale() const { return m_cmsScale; } - void setCMSScale(double cmsScale) { m_cmsScale=cmsScale; } - - double getDynamicScale() const { return m_dynamicScale; } - void setDynamicScale(double dynamicScale) { m_dynamicScale=dynamicScale; } - - - // set optimise flag on all sub grids - bool setOptimised(bool t=true) { - return m_optimised=t; - // for ( int iorder=0 ; iorder<2 ; iorder++ ) { - // for ( int iobs=0 ; iobs<Nobs() ; iobs++ ) m_grids[iorder][iobs]->setOptimised(t); - // } - } - - // find the number of words used for storage - int size() const; - - // get the cross sections - double& crossSection() { return m_total; } - double& crossSectionError() { return m_totalerror; } - - // double Lambda() const { return m_Lambda2; } - - // very lovely algebraic operators - grid& operator=(const grid& g); - grid& operator*=(const double& d); - grid& operator+=(const grid& g); - - /// test if grids have the same limits etc - bool operator==(const grid& g) const; - - // shouldn't have these, the grid is too large a structure - // to be passed in a return - // grid operator*(const double& d) const { return grid(*this)*=d; } - // grid operator+(const grid& g) const { return grid(*this)+=g; } - - void setDocumentation(const std::string& s); - void addDocumentation(const std::string& s); - - std::string getDocumentation() const { return m_documentation; } - std::string& getDocumentation() { return m_documentation; } - - - /// set the range of the observable bins, with an optional - /// scaling of the observable valuesfor channging units - void setBinRange(int ilower, int iupper, double xScaleFactor=1); - void setRange(double lower, double upper, double xScaleFactor=1); - - - /// add a correction as a std::vector - void addCorrection( std::vector<double>& v, const std::string& label="", bool combine=false ); - - - /// add a correction by histogram - void addCorrection(TH1D* h, const std::string& label="", double scale=1, bool combine=false ); - - - /// access the corrections - // const std::vector<std::vector<double> >& corrections() const { - const std::vector<correction>& corrections() const { - return m_corrections; - } - - /// get the correction labels - const std::vector<std::string >& correctionLabels() const { - return m_correctionLabels; - } - - /// will the corrections be applied? - bool getApplyCorrections() const { return m_applyCorrections; } - - bool setApplyCorrections(bool b) { - std::cout << "appl::grid bin-by-bin corrections will " - << ( b ? "" : "not " ) << "be applied" << std::endl; - return m_applyCorrections=b; - } - - /// apply corrections to a std::vector - void applyCorrections(std::vector<double>& v, std::vector<bool>& applied); - - - /// will a specific correction be applied? - bool getApplyCorrection(unsigned i) const { - if ( m_applyCorrections ) return true; - else if ( i<m_applyCorrection.size() ) return m_applyCorrection.at(i); - return false; - } - - bool setApplyCorrection(unsigned i, bool b) { - if ( i>=m_corrections.size() ) return false; - std::cout << "appl::grid bin-by-bin correction will " - << ( b ? "" : "not " ) << "be applied for correction " << i; - if ( m_correctionLabels[i]!="" ) std::cout << " (" << m_correctionLabels[i] << ")"; - std::cout << std::endl; - return m_applyCorrection[i]=b; - } - - /// apply corrections to a std::vector - bool applyCorrection(unsigned i, std::vector<double>& v); - - - /// set the ckm matrix values if need be - /// takes a 3x3 matrix with the format { { Vud, Vus, Vub }, { Vcd, Vcs, Vcb }, { Vtd, Vts, Vtb } } - void setckm( const std::vector<std::vector<double> >& ckm ); - - /// takes a flat 9 element vector (or c array) with the format { Vud, Vus, Vub, Vcd, Vcs, Vcb, Vtd, Vts, Vtb } - void setckm( const std::vector<double>& ckm ); - void setckm( const double* ckm ); - - - /// set the squared ckm matrix values if need be - /// the squared terms for eihter W+ or W- production - you probably should use setckm() - void setckm2( const std::vector<std::vector<double> >& ckm2 ); - - /// set the ckm matrix and squared ckm matrix values if need be - const std::vector<std::vector<double> >& getckm() const; - const std::vector<std::vector<double> >& getckm2() const; - - - /// flag custom convolution routines - - void sherpa() { m_type = SHERPA; std::cout << "appl::grid::sherpa() using SHERPA convolution" << std::endl; } - void amcatnlo() { m_type = AMCATNLO; std::cout << "appl::grid::amcatnlo() using aMC@NLO convolution" << std::endl; } - void standard() { m_type = STANDARD; std::cout << "appl::grid::standard() using standard convolution" << std::endl; } - - CALCULATION calculation() const { return m_type; } - - static std::string _calculation(CALCULATION C) { - switch (C) { - case STANDARD: - return "standard"; - case SHERPA: - return "sherpa"; - case AMCATNLO: - return "amcatnlo"; - case LAST_TYPE: - return "last_type"; // NB: shouldn't ever be used - } - return "unknown"; - } - - /// reduce number of subprocesses if possible - void shrink(const std::string& name, int ckmcharge=0); - - /// set bins to be combined after the convolution - void combine( std::vector<int>& v) { if ( (m_combine=v).size() ) combineReference(true); } - - /// set combine the be combined after the convolution - void combineReference(bool force=false); - - void combineBins(std::vector<double>& v, int power=1 ) const; - - double fx(double x) const; - double fy(double x) const; - - const appl_pdf* genpdf(int i) const { return m_genpdf[i]; } - - std::vector<double>& userdata() { return m_userdata; } - const std::vector<double>& userdata() const { return m_userdata; } - -protected: - - // internal common construct for the different types of constructor - void construct(int Nobs, - int NQ2=50, double Q2min=10000.0, double Q2max=25000000.0, int Q2order=4, - int Nx=50, double xmin=1e-5, double xmax=0.9, int xorder=3, - int order=2, - std::string transform="f" ); - -protected: - - /// std::string manipulators to parse the pdf names - - /// return chomped std::string - static std::string chomptoken(std::string& s1, const std::string& s2) - { - std::string s3 = ""; - std::string::size_type pos = s1.find(s2); - if ( pos != std::string::npos ) { - s3 = s1.substr(0, pos); - s1.erase(0, pos+1); - } - else { - s3 = s1.substr(0, s1.size()); - s1.erase(0, s1.size()+1); - } - return s3; - } - - static std::vector<std::string> parse(std::string s, const std::string& key) { - std::vector<std::string> clauses; - while ( s.size() ) clauses.push_back( chomptoken(s, key) ); - return clauses; - } - - /// get the required pdf combinations from those registered - void findgenpdf( std::string s ); - - /// add a generic pdf to the data base of registered pdfs - void addpdf( const std::string& s, const std::vector<int>& combinations=std::vector<int>() ); - - appl_pdf* genpdf(int i) { return m_genpdf[i]; } - -public: - - int subproc() const { return m_subproc; } - -protected: - - // histograms for saving the observable - TH1D* m_obs_bins; - TH1D* m_obs_bins_combined; - - // order in alpha_s of tree level contribution - int m_leading_order; - - // how many orders in the calculation, lo, nlo, nnlo etc - int m_order; - - // the actual weight grids themselves - igrid** m_grids[MAXGRIDS]; /// up to MAXGRIDS grids LO, NLO, NNLO, Real virtual, etc - - // total cross section qand uncertainty - double m_total; - double m_totalerror; - - // state variables - double m_run; - bool m_optimised; - bool m_trimmed; - - bool m_normalised; - - bool m_symmetrise; - - // transform and pdf combination tags - std::string m_transform; - std::string m_genpdfname; - - // pdf combination class - appl_pdf* m_genpdf[MAXGRIDS]; - - static const std::string m_version; - - double m_cmsScale; - - double m_dynamicScale; - - /// bin by bin correction factors - // std::vector<std::vector<double> > m_corrections; - std::vector<correction> m_corrections; - std::vector<std::string> m_correctionLabels; - - - /// should we apply the corrections? - bool m_applyCorrections; - - /// flag vector to determine whether each individual - /// correction should be applied - std::vector<bool> m_applyCorrection; - - std::string m_documentation; - - std::vector<double> m_ckmsum; - std::vector<std::vector<double> > m_ckm2; - std::vector<std::vector<double> > m_ckm; - - CALCULATION m_type; - - bool m_read; - - std::vector<int> m_combine; - - int m_subproc; - int m_bin; - - std::vector<double> m_userdata; - -}; - - -}; - -// shouldn't have this, grid is too large a structure -// grid operator*(const double& d, const appl::grid& g) { return g*d; } - -std::ostream& operator<<(std::ostream& s, const appl::grid& mygrid); - - - -#endif // __APPL_GRID_H diff --git a/reactions/APPLgrid/include/appl_grid/correction.h b/reactions/APPLgrid/include/appl_grid/correction.h deleted file mode 100644 index 647bbf86df3684e43fb86d9dec1ec9680efb2ad4..0000000000000000000000000000000000000000 --- a/reactions/APPLgrid/include/appl_grid/correction.h +++ /dev/null @@ -1,67 +0,0 @@ -// emacs: this is -*- c++ -*- -// -// @file correction.h -// class to store the multipliciative post processing -// corrections to be applied, only basic for the time -// but will be extended as appropriate -// -// Copyright (C) 2014 M.Sutton (sutt@cern.ch) -// -// $Id: correction.h, v0.0 Sun 23 Mar 2014 09:08:46 CET sutt $ - - -#ifndef CORRECTION_H -#define CORRECTION_H - -#include <iostream> -#include <vector> -#include <string> - - -// typedef std::vector<double> correction; - - -class correction { - -public: - - correction(const std::vector<double>& v, const std::string& s="" ) : mlabel(s), mv(v) { } - - virtual ~correction() { } - - std::string label() const { return mlabel; } - - unsigned size() const { return mv.size(); } - - double& operator[](int i) { return mv[i]; } - double operator[](int i) const { return mv[i]; } - - operator std::vector<double>&() { return mv; } - - correction operator=(const std::vector<double>& v) { mv=v; return *this; } - -private: - - std::string mlabel; - std::vector<double> mv; - -}; - - -// inline std::ostream& operator<<( std::ostream& s, const correction& /* _c */ ) { -// return s; -// } - - - -#endif // CORRECTION_H - - - - - - - - - - diff --git a/reactions/APPLgrid/src/ReactionAPPLgrid.cc b/reactions/APPLgrid/src/ReactionAPPLgrid.cc index 99396b5141c74c55b29ad91879c46f9e61776c98..a50fa91704f9b922823c798df3bf66901a244634 100644 --- a/reactions/APPLgrid/src/ReactionAPPLgrid.cc +++ b/reactions/APPLgrid/src/ReactionAPPLgrid.cc @@ -29,7 +29,6 @@ void ReactionAPPLgrid::setDatasetParameters(int dataSetID, map<string,string> pa std::string token; while(std::getline(ss, token, ',')) { - //std::cout << token << '\n'; std::shared_ptr<appl::grid> g(new appl::grid(token)); g->trim(); _grids[dataSetID].push_back(g); diff --git a/tools/draw/include/PdfData.h b/tools/draw/include/PdfData.h index 5df80031e728d1f8539ab6caefd6b1a281086fed..0e053bd9feab60ceeba9d4c32d2a18c852397b21 100644 --- a/tools/draw/include/PdfData.h +++ b/tools/draw/include/PdfData.h @@ -11,6 +11,7 @@ using namespace std; enum pdferr {None, AsymHess, SymHess, MC}; enum pdftype{uv=0, dv, g, Sea, ubar, dbar, s, Rs, c, b, dbarminubar, uvmindv, U, D, Ubar, Dbar, goversea, doveru, dbaroverubar, dvoveruv, rs, photon, SeaOverGlue, photonOverGlue}; +//enum pdftype{uv=0, dv, g, Sea, ubar, dbar, s, Rs, c, b, dbarminubar, uvmindv, U, D, Ubar, Dbar, goversea, doveru, dbaroverubar, dvoveruv, rs, photon, SeaOverGlue, photonOverGlue, uvplusdv, uvplusdvplusSea,uvplusdvminSea}; extern vector <pdftype> pdfs; extern vector <string> pdflabels; diff --git a/tools/draw/src/PdfData.cc b/tools/draw/src/PdfData.cc index 034a3ebe5847cb2c94b0d423c686448d56980ecd..526d13b9bf8e83ad01661c385bb1d0a58df1824c 100644 --- a/tools/draw/src/PdfData.cc +++ b/tools/draw/src/PdfData.cc @@ -20,12 +20,15 @@ #include "FileOpener.h" //pdf type -pdftype pdfts[] = {uv, dv, g, Sea, ubar, dbar, s, Rs, c, b, dbarminubar, uvmindv, U, D, Ubar, Dbar, goversea, doveru, dbaroverubar, dvoveruv,rs,photon,SeaOverGlue, photonOverGlue }; +pdftype pdfts[] = {uv, dv, g, Sea, ubar, dbar, s, Rs, c, b, dbarminubar, uvmindv, U, D, Ubar, Dbar, goversea, doveru, dbaroverubar, dvoveruv,rs,photon,SeaOverGlue, photonOverGlue}; +//pdftype pdfts[] = {uv, dv, g, Sea, ubar, dbar, s, Rs, c, b, dbarminubar, uvmindv, U, D, Ubar, Dbar, goversea, doveru, dbaroverubar, dvoveruv,rs,photon,SeaOverGlue, photonOverGlue, uvplusdv, uvplusdvplusSea,uvplusdvminSea }; //pdf labels string pdflab[] = {"u_{V}", "d_{V}", "g", "#Sigma", "#bar{u}", "#bar{d}", "s", "(s+#bar{s})/(#bar{u}+#bar{d})", "c", "b", "#bar{d}-#bar{u}", "d_{V}-u_{V}", "U", "D", "#bar{U}", "#bar{D}", "g/#Sigma", "d/u", "#bar{d}/#bar{u}", "d_{V}/u_{V}","rs","#gamma","#Sigma/g","#gamma/g"}; +// "d/u", "#bar{d}/#bar{u}", "d_{V}/u_{V}","rs","#gamma","#Sigma/g","#gamma/g","uv+dv","uv+dv+2#Sigma","uv+dv-2#Sigma"}; //pdf filenames string pdffil[] = {"uv", "dv", "g", "Sea", "ubar", "dbar", "s", "Rs", "c", "b", "dbar-ubar", "uv-dv", "U", "D", "UBar", "DBar", "goversea", "doveru", "dbaroverubar", "dvoveruv","rs","ph","sg","gg" +//string pdffil[] = {"uv", "dv", "g", "Sea", "ubar", "dbar", "s", "Rs", "c", "b", "dbar-ubar", "uv-dv", "U", "D", "UBar", "DBar", "goversea", "doveru", "dbaroverubar", //"dvoveruv","rs","ph","sg","gg","uv+dv","uv+dv+2Sea","uv+dv-2Sea" }; vector <pdftype> pdfs(pdfts, pdfts + sizeof(pdfts) / sizeof(pdftype)); @@ -154,6 +157,27 @@ Pdf::Pdf(string filename) : Q2value(0), NxValues(0), NPdfs(0), Xmin(0), Xmax(0) else tablemap[photonOverGlue].push_back(0); + /*PdfTypes.push_back(uvplusdv); NPdfs++; + for (int ix = 0; ix < NxValues; ix++) + if (tablemap[g][ix] != 0) + tablemap[uvplusdv].push_back(tablemap[dv][ix]+tablemap[uv][ix]); + else + tablemap[uvplusdv].push_back(0); + + PdfTypes.push_back(uvplusdvplusSea); NPdfs++; + for (int ix = 0; ix < NxValues; ix++) + if (tablemap[g][ix] != 0) + tablemap[uvplusdvplusSea].push_back(tablemap[dv][ix]+tablemap[uv][ix]+2.0*tablemap[Sea][ix]); + else + tablemap[uvplusdvplusSea].push_back(0); + + PdfTypes.push_back(uvplusdvminSea); NPdfs++; + for (int ix = 0; ix < NxValues; ix++) + if (tablemap[g][ix] != 0) + tablemap[uvplusdvminSea].push_back(tablemap[dv][ix]+tablemap[uv][ix]-2.0*tablemap[Sea][ix]); + else + tablemap[uvplusdvminSea].push_back(0);*/ + }