Update box matrix normalization and decouple K matrix model from BoxQuantization

README

@@ -19,6 +19,8 @@ An overview of this software is given in the paper below:
       doi            = "10.1016/j.nuclphysb.2017.09.014"
 }
 
+**NOTE:** The normalization of box matrix elements was changed in https://github.com/ebatz/TwoHadronsInBox/commit/f5f372f8ffc09f217cdb2796685b1e8551251be5 compared to the above reference. The consequences for the K matrix normalization are explained in the commit message.
+
 
 See comments in the header files for information on each class
 and required XML forms.

build/Makefile

@@ -36,7 +36,8 @@ CXXFLAGS = -O3 -std=c++11 -Wall -I../source
 DFLAGS = -DAR_LMAX=$(AR_LMAX_VALUE) -DOA_LMAX=$(OA_LMAX_VALUE) \
          -DPD_LMAX=$(PD_LMAX_VALUE) -DCD_LMAX=$(CD_LMAX_VALUE) \
          -DAR_SX2MAX=$(AR_SX2MAX_VALUE) -DOA_SX2MAX=$(OA_SX2MAX_VALUE) \
-         -DPD_SX2MAX=$(PD_SX2MAX_VALUE) -DCD_SX2MAX=$(CD_SX2MAX_VALUE)
+         -DPD_SX2MAX=$(PD_SX2MAX_VALUE) -DCD_SX2MAX=$(CD_SX2MAX_VALUE) \
+	 -fPIC
 LDFLAGS =
 LIBS=-llapack
 
@@ -63,6 +64,7 @@ INCS = \
   box_quant.h      \
   cmframe.h        \
   fit_forms.h      \
+  K_matrix_base.h  \
   K_matrix_calc.h  \
   K_matrix_info.h  \
   matrix.h         \
@@ -76,6 +78,9 @@ OBJS = $(SRCS:.cc=.o)
 
 all:  driver1 driver2
 
+lib: $(OBJS)
+	$(CXX) $(CXXFLAGS) -shared $(LDFLAGS) -o libBox.so $(OBJS) $(LDFLAGS)
+
 driver%: driver%.cc $(OBJS)
 	$(CXX) $(CXXFLAGS) $(LDFLAGS) -o driver$* ../source/driver$*.cc $(OBJS) $(LIBS) 
 

source/K_matrix_base.h

@@ -0,0 +1,25 @@
+#ifndef K_MATRIX_BASE_H
+#define K_MATRIX_BASE_H
+
+#include <map>
+
+
+class KtildeMatrixBase
+{
+
+ public:
+
+    virtual ~KtildeMatrixBase() {};
+
+    virtual double calculate(uint Jtimestwo, 
+                     uint Lp, uint Sptimestwo, uint chanp,
+                     uint L, uint Stimestwo, uint chan,
+                     double Ecm_over_mref) const = 0;
+
+    virtual bool isZero(uint Jtimestwo, 
+                uint Lp, uint Sptimestwo, uint chanp,
+                uint L, uint Stimestwo, uint chan) const = 0;
+
+};
+
+#endif

source/K_matrix_calc.cc

@@ -172,18 +172,25 @@ KtildeMatrixCalculator::KtildeMatrixCalculator(const std::list<std::pair<KElemen
 }
 
 
+void KtildeMatrixCalculator::setKtildeParameters(const std::vector<double>& kappa_params)
+{
+ if (kappa_params.size()!=getNumberOfParameters())
+    throw(std::invalid_argument("Could not set KtildeParameters"));
+ m_kappa_params=kappa_params;
+}
+
+
 double KtildeMatrixCalculator::calculate(uint Jtimestwo, 
                        uint Lp, uint Sptimestwo, uint chanp,
                        uint L, uint Stimestwo, uint chan, 
-                       const vector<double>& kappa,
                        double Ecm_over_mref) const
 {
- if (kappa.size()<m_paraminfo.size())
+ if (m_kappa_params.size()<m_paraminfo.size())
     throw(std::invalid_argument("Insufficient number of parameter values in KtildeMatrixCalculator::calculate"));
  KElementInfo kelem(Jtimestwo,Lp,Sptimestwo,chanp,L,Stimestwo,chan);
  std::map<KElementInfo,FitForm*>::const_iterator it=m_fit.find(kelem);
  if (it==m_fit.end()) return 0.0;
- return (it->second)->evaluate(kappa,Ecm_over_mref);
+ return (it->second)->evaluate(m_kappa_params,Ecm_over_mref);
 }
 
 
@@ -343,18 +350,25 @@ KtildeInverseCalculator::KtildeInverseCalculator(const std::list<std::pair<KElem
 }
 
 
+void KtildeInverseCalculator::setKtildeParameters(const std::vector<double>& kappa_params)
+{
+ if (kappa_params.size()!=getNumberOfParameters())
+    throw(std::invalid_argument("Could not set KtildeInvParameters"));
+ m_kappa_params=kappa_params;
+}
+
+
 double KtildeInverseCalculator::calculate(uint Jtimestwo, 
                        uint Lp, uint Sptimestwo, uint chanp,
                        uint L, uint Stimestwo, uint chan, 
-                       const vector<double>& kappa,
                        double Ecm_over_mref) const
 {
- if (kappa.size()<m_paraminfo.size())
+ if (m_kappa_params.size()<m_paraminfo.size())
     throw(std::invalid_argument("Insufficient number of parameter values in KtildeInverseCalculator::calculate"));
  KElementInfo kelem(Jtimestwo,Lp,Sptimestwo,chanp,L,Stimestwo,chan);
  std::map<KElementInfo,FitForm*>::const_iterator it=m_fit.find(kelem);
  if (it==m_fit.end()) return 0.0;
- return (it->second)->evaluate(kappa,Ecm_over_mref);
+ return (it->second)->evaluate(m_kappa_params,Ecm_over_mref);
 }
 
 

source/K_matrix_calc.h

@@ -4,6 +4,7 @@
 #include "xml_handler.h"
 #include "fit_forms.h"
 #include "K_matrix_info.h"
+#include "K_matrix_base.h"
 #include <map>
 
 // ***********************************************************************
@@ -70,12 +71,13 @@
 
 
 
-class KtildeMatrixCalculator
+class KtildeMatrixCalculator : public KtildeMatrixBase
 {
 
     std::map<KElementInfo,FitForm*> m_fit;
     std::vector<KFitParamInfo> m_paraminfo;
     std::map<KFitParamInfo,uint> m_paramindices;
+    std::vector<double> m_kappa_params;
 
         // prevent copying, no default
 
@@ -109,10 +111,11 @@ class KtildeMatrixCalculator
     void output(XMLHandler& xmlout) const;  // XML output
 
 
+    void setKtildeParameters(const std::vector<double>& kappa_params);
+
     double calculate(uint Jtimestwo, 
                      uint Lp, uint Sptimestwo, uint chanp,
                      uint L, uint Stimestwo, uint chan,
-                     const std::vector<double>& kappa_params,
                      double Ecm_over_mref) const;
 
     bool isZero(uint Jtimestwo, 
@@ -136,12 +139,13 @@ class KtildeMatrixCalculator
 
 
 
-class KtildeInverseCalculator
+class KtildeInverseCalculator : public KtildeMatrixBase
 {
 
     std::map<KElementInfo,FitForm*> m_fit;
     std::vector<KFitParamInfo> m_paraminfo;
     std::map<KFitParamInfo,uint> m_paramindices;
+    std::vector<double> m_kappa_params;
 
          // disallow copying, no default
 
@@ -174,10 +178,11 @@ class KtildeInverseCalculator
     void output(XMLHandler& xmlout) const;  // XML output
 
 
+    void setKtildeParameters(const std::vector<double>& kappa_params);
+
     double calculate(uint Jtimestwo, 
                      uint Lp, uint Sptimestwo, uint chanp,
                      uint L, uint Stimestwo, uint chan,
-                     const std::vector<double>& kappa_params,
                      double Ecm_over_mref) const;
 
     bool isZero(uint Jtimestwo,