modified: ../core/base/include/Physics.h

	modified:   ../core/base/src/Physics.cxx
	modified:   asymmetries/include/VCSAsFromPDFs.h
	modified:   structurefunctions/include/CompositeSFs.h
	modified:   structurefunctions/include/SFsFromPDFs.h
	modified:   structurefunctions/include/SpinStructureFunctions.h
	modified:   structurefunctions/src/SpinStructureFunctions.cxx

src/core/base/include/Physics.h

@@ -38,6 +38,24 @@ namespace insane {
 
     /** @}*/
 
+    /** \addtogroup Deep Inelastic Scattering 
+     *  @{
+     */
+
+    /** R \f$ = \sigma_L/\sigma_T\f$ (ratio of L/T cross sections).
+     *
+     * Double_t res = F2p(x, Q2) / (2.0 * x * F1p(x, Q2)) ;
+     * res = res * (1 + 4.0 * M_p/GeV * M_p/GeV * x * x / Q2) - 1.0;
+     *
+     */
+    double R(double F1, double F2, double x, double Q2);
+
+    /** Soffer-Teryaev  Positivity bound.
+     *
+     *  \f$ |A_2| < \sqrt{R(1+A_1)/2} \f$
+     */
+    double A2_PositivityBound(double F_1, double F_2, double A_1, double x, double Q2);
+    double A2_PositivityBound(double R, double A_1, double x, double Q2);
 
     /** Wandzura-Wilczek function form.
      *
@@ -46,6 +64,9 @@ namespace insane {
      */
     double WW(std::function<double(double)> func, double x, double max=1.0);
 
+    /** @}*/
+
+
     namespace TMCs {
       // Target Mass correction function
       // J. Phys. G 35, 053101 (2008)

src/core/base/src/Physics.cxx

@@ -22,6 +22,26 @@ namespace insane {
       return (GM - GE)/(1.0+t);
     }
 
+    double R(double F1, double F2, double x, double Q2)
+    {
+      double res = F2 / (2.0 * x * F1) ;
+      res = res * (1 + 4.0 * (M_p/GeV)*(M_p/GeV)* x * x / Q2) - 1.0;
+      return res;
+    }
+
+    // Soffer-Teryaev  Positivity bound.
+    // |A_2| < \sqrt{R(1+A_1)/2}
+    double A2_PositivityBound(double F_1, double F_2, double A_1, double x, double Q2){
+      //double F_1 = F1(updfs, target, x, Q2);
+      //double F_2 = F2(updfs, target, x, Q2);
+      double RLT = insane::physics::R(F_1,F_2,x,Q2);
+      //double A_2 = A2_NoTMC(updfs,ppdfs,target,x,Q2);
+      return TMath::Sqrt(RLT*(1.0+A_1)/2.0 );
+    }
+    double A2_PositivityBound(double R, double A_1, double x, double Q2){
+      return TMath::Sqrt( R*(1.0+A_1)/2.0 );
+    }
+
     double WW(std::function<double(double)> func, double x, double max)
     {
       return( -1.0*func(x) + insane::integrate::simple(func,x,max) );

src/physics/asymmetries/include/VCSAsFromPDFs.h

@@ -3,6 +3,7 @@
 
 #include "VirtualComptonScatteringAsymmetries.h"
 #include "SSFsFromPDFs.h"
+#include "SFsFromPDFs.h"
 #include "TMath.h"
 #include "InSANEMathFunc.h"
 
@@ -104,6 +105,16 @@ namespace insane {
       return( gam*(g_1 + g_2)/F_1 );
     } 
 
+    // Soffer-Teryaev  Positivity bound.
+    // |A_2| < \sqrt{R(1+A_1)/2}
+    template< class...T1, class...T2 >
+    double A2_PositivityBound(const std::tuple<T1...>& updfs, const std::tuple<T2...>& ppdfs, Nuclei target, double x, double Q2){
+      double F_1 = F1(updfs, target, x, Q2);
+      double F_2 = F2(updfs, target, x, Q2);
+      double RLT = insane::physics::R(F_1,F_2,x,Q2);
+      double A_2 = A2_NoTMC(updfs,ppdfs,target,x,Q2);
+      return TMath::Sqrt( RLT*(1.0+A_2)/2.0 );
+    }
 
     template<class...> 
       class VCSAsFromPDFs{};

src/physics/structurefunctions/include/CompositeSFs.h

@@ -5,7 +5,7 @@
 #include "Physics.h"
 #include "NMC95_SFs.h"
 #include "F1F209_SFs.h"
-#include "MAID_VPACs.h"
+//#include "MAID_VPACs.h"
 #include "SFsFromVPACs.h"
 #include "SFsFromPDFs.h"
 #include "Stat2015_UPDFs.h"
@@ -13,7 +13,7 @@
 namespace insane {
   namespace physics {
 
-    class  CompositeSFs : public StructureFunctions { 
+    class  CompositeSFs : public StructureFunctions2 { 
       protected:
         NMC95_SFs   fAsf;
         //SFsFromVPACs<MAID_VPACs> fBsf;

src/physics/structurefunctions/include/SFsFromPDFs.h

@@ -63,6 +63,15 @@ namespace insane {
     } 
     //______________________________________________________________________________
 
+    //template<class T>
+    //double R(const T& pdfs, Nuclei target, double x, double Q2){
+    //  // R=sigma_L/sigma_T
+    //  double res = pdfs.F2p(x,Q2)/(2.0*x*pdfs.F1p(x, Q2));
+    //  res = res * (1 + 4.0*(M_p/GeV)*(M_p/GeV)*x*x/Q2)-1.0;
+    //  if(TMath::IsNaN(res) ) return 0.0;
+    //  return res;
+    //}
+
     /** StructureFunction implementation for polarized PDFs. 
      *  Although we don't need this class to calculate the SSFs, it provides 
      *  some usefulness when using with other structure functions.

src/physics/structurefunctions/include/SpinStructureFunctions.h

@@ -97,7 +97,7 @@ namespace insane {
         //virtual double d2p_Twist3_TMC(double Q2,double x1 = 0.01,double x2 = 0.99) ;
 
         virtual double d2p_tilde(double Q2,double x1 = 0.01,double x2 = 0.99) ;
-        //virtual double d2n_tilde(double Q2,double x1 = 0.01,double x2 = 0.99) ;
+        virtual double d2n_tilde(double Q2,double x1 = 0.01,double x2 = 0.99) ;
 
         virtual double d2p_tilde_TMC(double Q2,double x1 = 0.01,double x2 = 0.99) ;
 

src/physics/structurefunctions/src/SpinStructureFunctions.cxx

@@ -252,6 +252,23 @@ namespace insane {
     }
     //_____________________________________________________________________________
 
+    double SpinStructureFunctions::d2n_tilde(double Q2,double x1,double x2)
+    {
+      // d2 tilde (no elastic contribution) 
+      //double x1     = 0.1;
+      //double x2     = 0.9;
+      double result = 0.0;
+      Int_t N         = 100;
+      double dx     = (x2 - x1) / ((double)N);
+      // quick simple integration
+      for (int i = 0; i < N; i++) {
+        double x = x1 + dx*double(i);
+        result += (dx*(x*x)*( 2.0*g1n(x,Q2) + 3.0*g2n(x,Q2) ));
+      }
+      return(result);
+    }
+    //_____________________________________________________________________________
+
     double SpinStructureFunctions::d2p_tilde_TMC(double Q2,double x1,double x2)
     {
       // d2 tilde (no elastic contribution)