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);*/
+
 
 }