From 35dff2fc3fd0397b7922793c6ef63d5f0d4b868e Mon Sep 17 00:00:00 2001 From: Ben Horz Date: Wed, 18 Apr 2018 16:04:11 +0200 Subject: [PATCH 1/4] Add build option for shared library. --- build/Makefile | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/build/Makefile b/build/Makefile index 339db6d..d7b334b 100644 --- a/build/Makefile +++ b/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 @@ -76,6 +77,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) From de71b5791b890c0a6ad060ddab66cf8a8614bc9f Mon Sep 17 00:00:00 2001 From: Ben Horz Date: Wed, 18 Apr 2018 16:13:44 +0200 Subject: [PATCH 2/4] Decouple K matrix parametrization from BoxQuantization. Previously the K matrix model was intertwined with the BoxQuantizationclass, e.g. fit parameters were accessed through BoxQuantization. Instead, the BoxQuantization class should not care about how a K Matrix model is implemented as long as it provides the two functions calculate(..) isZero(..) now defined pure virtual in the new KtildeMatrixBase base class. A K-matrix model is implemented by subclassing KtildeMatrixBase and providing these functions. Model parameters are handled exclusively by the K matrix model provided by the user. A BoxQuantization object is now constructed by passing in an object derived from KtildeMatrixBase and a boolean indicating whether this is a model for Ktilde or KtildeInverse. The example driver routines have been updated to reflect the changes in interface. --- build/Makefile | 1 + source/K_matrix_base.h | 25 ++++++ source/K_matrix_calc.cc | 26 ++++-- source/K_matrix_calc.h | 13 ++- source/box_quant.cc | 181 +++++++++------------------------------- source/box_quant.h | 47 ++--------- source/driver1.cc | 27 ++++-- source/driver2.cc | 6 +- 8 files changed, 126 insertions(+), 200 deletions(-) create mode 100644 source/K_matrix_base.h diff --git a/build/Makefile b/build/Makefile index d7b334b..723d5d6 100644 --- a/build/Makefile +++ b/build/Makefile @@ -64,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 \ diff --git a/source/K_matrix_base.h b/source/K_matrix_base.h new file mode 100644 index 0000000..ab5ce7d --- /dev/null +++ b/source/K_matrix_base.h @@ -0,0 +1,25 @@ +#ifndef K_MATRIX_BASE_H +#define K_MATRIX_BASE_H + +#include + + +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 diff --git a/source/K_matrix_calc.cc b/source/K_matrix_calc.cc index 754b26c..d2e20f9 100644 --- a/source/K_matrix_calc.cc +++ b/source/K_matrix_calc.cc @@ -172,18 +172,25 @@ KtildeMatrixCalculator::KtildeMatrixCalculator(const std::list& 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& kappa, double Ecm_over_mref) const { - if (kappa.size()::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& 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& kappa, double Ecm_over_mref) const { - if (kappa.size()::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); } diff --git a/source/K_matrix_calc.h b/source/K_matrix_calc.h index 8bc6fd5..cf5ea60 100644 --- a/source/K_matrix_calc.h +++ b/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 // *********************************************************************** @@ -70,12 +71,13 @@ -class KtildeMatrixCalculator +class KtildeMatrixCalculator : public KtildeMatrixBase { std::map m_fit; std::vector m_paraminfo; std::map m_paramindices; + std::vector 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& kappa_params); + double calculate(uint Jtimestwo, uint Lp, uint Sptimestwo, uint chanp, uint L, uint Stimestwo, uint chan, - const std::vector& 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 m_fit; std::vector m_paraminfo; std::map m_paramindices; + std::vector 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& kappa_params); + double calculate(uint Jtimestwo, uint Lp, uint Sptimestwo, uint chanp, uint L, uint Stimestwo, uint chan, - const std::vector& kappa_params, double Ecm_over_mref) const; bool isZero(uint Jtimestwo, diff --git a/source/box_quant.cc b/source/box_quant.cc index d40faa4..e70fa59 100644 --- a/source/box_quant.cc +++ b/source/box_quant.cc @@ -200,33 +200,14 @@ void BoxQuantBasisState::output(XMLHandler& xmlout) const // * BoxQuantization::BoxQuantization(XMLHandler& xmlin, - KtildeMatrixCalculator *Kmatptr) - : m_Kmat(Kmatptr), m_Kinv(0) + KtildeMatrixBase *Kmatptr, + bool KInvMode) + : m_Kptr(Kmatptr), m_KInvMode(KInvMode) { xmlinitialize(xmlin); } -BoxQuantization::BoxQuantization(XMLHandler& xmlin, - KtildeInverseCalculator *Kinvptr) - : m_Kmat(0), m_Kinv(Kinvptr) -{ - xmlinitialize(xmlin); -} - - -BoxQuantization::BoxQuantization(XMLHandler& xmlin, - KtildeMatrixCalculator *Kmatptr, - KtildeInverseCalculator *Kinvptr) - : m_Kmat(Kmatptr), m_Kinv(Kinvptr) -{ - if ((m_Kmat!=0)&&(m_Kinv!=0)) - throw(std::invalid_argument(string("One of the pointers must be null in BoxQuantization constructor with ") - +string("KtildeMatrixCalculator and KtildeInverseCalculator pointers"))); - xmlinitialize(xmlin); -} - - void BoxQuantization::xmlinitialize(XMLHandler& xmlin) { try{ @@ -254,24 +235,10 @@ BoxQuantization::BoxQuantization(const std::string& mom_ray, uint mom_int_sq, const std::string& lgirrep, const std::vector& chan_infos, const std::vector lmaxes, - KtildeMatrixCalculator *Kmatptr) - : m_lgirrep(lgirrep), m_momray(mom_ray), m_decay_infos(chan_infos), - m_Lmaxes(lmaxes), m_Kmat(Kmatptr), m_Kinv(0) -{ - set_dvector(mom_int_sq); - if (m_decay_infos.size()!=m_Lmaxes.size()) - throw(std::invalid_argument("Size mismatch between decay infos and L-maxes")); - initialize(); -} - - -BoxQuantization::BoxQuantization(const std::string& mom_ray, uint mom_int_sq, - const std::string& lgirrep, - const std::vector& chan_infos, - const std::vector lmaxes, - KtildeInverseCalculator *Kinvptr) + KtildeMatrixBase *Kmatptr, + bool KInvMode) : m_lgirrep(lgirrep), m_momray(mom_ray), m_decay_infos(chan_infos), - m_Lmaxes(lmaxes), m_Kmat(0), m_Kinv(Kinvptr) + m_Lmaxes(lmaxes), m_Kptr(Kmatptr), m_KInvMode(KInvMode) { set_dvector(mom_int_sq); if (m_decay_infos.size()!=m_Lmaxes.size()) @@ -364,8 +331,6 @@ void BoxQuantization::initialize() m_masses1[k]=m_masses2[k]=1.0; // arbitrary default values // set up the basis states setup_basis(); - // initialize size of parameters - m_kappa_params.resize(getNumberOfKtildeParameters()); } @@ -402,12 +367,6 @@ void BoxQuantization::setMassesOverRef(uint channel_index, double mass1_over_ref } -void BoxQuantization::setKtildeParameters(std::vector kappa_params) -{ - if (kappa_params.size()!=getNumberOfKtildeParameters()) - throw(std::invalid_argument("Could not set KtildeParameters")); - m_kappa_params=kappa_params; -} string BoxQuantization::output(int indent) const @@ -460,55 +419,6 @@ string BoxQuantization::outputBasis(int indent) const } -void BoxQuantization::outputKFitParams(XMLHandler& xmlout) const // XML output -{ - xmlout.set_root("BoxQuantizationKFitParams"); - const std::vector& fref=getFitParameterInfos(); - for (std::vector::const_iterator it=fref.begin();it!=fref.end();it++){ - XMLHandler xmlk; - it->output(xmlk); - xmlout.put_child(xmlk);} -} - - -std::string BoxQuantization::outputKFitParams(int indent) const // XML output -{ - XMLHandler xmlout; - outputKFitParams(xmlout); - return xmlout.output(indent); -} - - - - - -const std::vector& BoxQuantization::getFitParameterInfos() const -{ - if (m_Kmat) return m_Kmat->getFitParameterInfos(); - else return m_Kinv->getFitParameterInfos(); -} - - -int BoxQuantization::getParameterIndex(const KFitParamInfo& kinfo) const // returns -1 if not found -{ - if (m_Kmat) return m_Kmat->getParameterIndex(kinfo); - else return m_Kinv->getParameterIndex(kinfo); -} - - -double BoxQuantization::getParameterValue(const KFitParamInfo& kinfo) const -{ - if (m_Kmat) return m_kappa_params.at(m_Kmat->getParameterIndex(kinfo)); - else return m_kappa_params.at(m_Kinv->getParameterIndex(kinfo)); -} - -set BoxQuantization::getElementInfos() const -{ - if (m_Kmat) return m_Kmat->getElementInfos(); - else return m_Kinv->getElementInfos(); -} - - string BoxQuantization::getLGIrrepParityFlip() { if (m_momray=="ar"){ @@ -584,8 +494,7 @@ void BoxQuantization::setup_basis() bqm.getColumnJtimestwo(),bqm.getColumnL(),bqm.getColumnNocc()));}}} // look for states to exclude due to K-matrix set exclusions; - if (m_Kmat!=0) exclusions=find_excluded_states_from_ktilde(m_Kmat); - else exclusions=find_excluded_states_from_ktilde(m_Kinv); + exclusions=find_excluded_states_from_ktilde(m_Kptr); // if ALL states excluded, this means Kmatrix is zero; end user needs to know // which K matrix elements should be set; throw this information if (exclusions.size()==m_basis.size()){ @@ -649,105 +558,93 @@ void BoxQuantization::getBoxMatrixFromEcm(double Ecm_over_mref, CMatrix& B) void BoxQuantization::getKtildeFromElab(double Elab_over_mref, RealSymmetricMatrix& Ktilde) { - if (m_Kmat==0) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); + if (m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); RMatrix dummy; - get_ktilde_matrix(Elab_over_mref,Ktilde,dummy,true,true,m_Kmat); + get_ktilde_matrix(Elab_over_mref,Ktilde,dummy,true,true,m_Kptr); } void BoxQuantization::getKtildeFromElab(double Elab_over_mref, RMatrix& Ktilde) { - if (m_Kmat==0) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); + if (m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); RealSymmetricMatrix dummy; - get_ktilde_matrix(Elab_over_mref,dummy,Ktilde,true,false,m_Kmat); + get_ktilde_matrix(Elab_over_mref,dummy,Ktilde,true,false,m_Kptr); } void BoxQuantization::getKtildeFromEcm(double Ecm_over_mref, RealSymmetricMatrix& Ktilde) { - if (m_Kmat==0) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); + if (m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); RMatrix dummy; - get_ktilde_matrix(Ecm_over_mref,Ktilde,dummy,false,true,m_Kmat); + get_ktilde_matrix(Ecm_over_mref,Ktilde,dummy,false,true,m_Kptr); } void BoxQuantization::getKtildeFromEcm(double Ecm_over_mref, RMatrix& Ktilde) { - if (m_Kmat==0) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); + if (m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktilde in Ktildeinverse mode")); RealSymmetricMatrix dummy; - get_ktilde_matrix(Ecm_over_mref,dummy,Ktilde,false,false,m_Kmat); + get_ktilde_matrix(Ecm_over_mref,dummy,Ktilde,false,false,m_Kptr); } void BoxQuantization::getKtildeinvFromElab(double Elab_over_mref, RealSymmetricMatrix& Ktildeinv) { - if (m_Kinv==0) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); + if (!m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); RMatrix dummy; - get_ktilde_matrix(Elab_over_mref,Ktildeinv,dummy,true,true,m_Kinv); + get_ktilde_matrix(Elab_over_mref,Ktildeinv,dummy,true,true,m_Kptr); } void BoxQuantization::getKtildeinvFromElab(double Elab_over_mref, RMatrix& Ktildeinv) { - if (m_Kinv==0) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); + if (!m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); RealSymmetricMatrix dummy; - get_ktilde_matrix(Elab_over_mref,dummy,Ktildeinv,true,false,m_Kinv); + get_ktilde_matrix(Elab_over_mref,dummy,Ktildeinv,true,false,m_Kptr); } void BoxQuantization::getKtildeinvFromEcm(double Ecm_over_mref, RealSymmetricMatrix& Ktildeinv) { - if (m_Kinv==0) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); + if (!m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); RMatrix dummy; - get_ktilde_matrix(Ecm_over_mref,Ktildeinv,dummy,false,true,m_Kinv); + get_ktilde_matrix(Ecm_over_mref,Ktildeinv,dummy,false,true,m_Kptr); } void BoxQuantization::getKtildeinvFromEcm(double Ecm_over_mref, RMatrix& Ktildeinv) { - if (m_Kinv==0) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); + if (!m_KInvMode) throw(std::runtime_error("Cannot evaluate Ktildeinverse in Ktilde mode")); RealSymmetricMatrix dummy; - get_ktilde_matrix(Ecm_over_mref,dummy,Ktildeinv,false,false,m_Kinv); + get_ktilde_matrix(Ecm_over_mref,dummy,Ktildeinv,false,false,m_Kptr); } void BoxQuantization::getKtildeOrInverseFromElab(double Elab_over_mref, RealSymmetricMatrix& KtildeOrInverse) { RMatrix dummy; - if (m_Kmat!=0) - get_ktilde_matrix(Elab_over_mref,KtildeOrInverse,dummy,true,true,m_Kmat); - else - get_ktilde_matrix(Elab_over_mref,KtildeOrInverse,dummy,true,true,m_Kinv); + get_ktilde_matrix(Elab_over_mref,KtildeOrInverse,dummy,true,true,m_Kptr); } void BoxQuantization::getKtildeOrInverseFromElab(double Elab_over_mref, RMatrix& KtildeOrInverse) { RealSymmetricMatrix dummy; - if (m_Kmat!=0) - get_ktilde_matrix(Elab_over_mref,dummy,KtildeOrInverse,true,false,m_Kmat); - else - get_ktilde_matrix(Elab_over_mref,dummy,KtildeOrInverse,true,false,m_Kinv); + get_ktilde_matrix(Elab_over_mref,dummy,KtildeOrInverse,true,false,m_Kptr); } void BoxQuantization::getKtildeOrInverseFromEcm(double Ecm_over_mref, RealSymmetricMatrix& KtildeOrInverse) { RMatrix dummy; - if (m_Kmat!=0) - get_ktilde_matrix(Ecm_over_mref,KtildeOrInverse,dummy,false,true,m_Kmat); - else - get_ktilde_matrix(Ecm_over_mref,KtildeOrInverse,dummy,false,true,m_Kinv); + get_ktilde_matrix(Ecm_over_mref,KtildeOrInverse,dummy,false,true,m_Kptr); } void BoxQuantization::getKtildeOrInverseFromEcm(double Ecm_over_mref, RMatrix& KtildeOrInverse) { RealSymmetricMatrix dummy; - if (m_Kmat!=0) - get_ktilde_matrix(Ecm_over_mref,dummy,KtildeOrInverse,false,false,m_Kmat); - else - get_ktilde_matrix(Ecm_over_mref,dummy,KtildeOrInverse,false,false,m_Kinv); + get_ktilde_matrix(Ecm_over_mref,dummy,KtildeOrInverse,false,false,m_Kptr); } @@ -880,16 +777,16 @@ void BoxQuantization::assign_matrices(double E_over_mref, bool Elab, ComplexHerm { if (Elab){ getBoxMatrixFromElab(E_over_mref,B); - if (m_Kmat!=0) - getKtildeFromElab(E_over_mref,Kv); + if (m_KInvMode) + getKtildeinvFromElab(E_over_mref,Kv); else - getKtildeinvFromElab(E_over_mref,Kv);} + getKtildeFromElab(E_over_mref,Kv);} else{ getBoxMatrixFromEcm(E_over_mref,B); - if (m_Kmat!=0) - getKtildeFromEcm(E_over_mref,Kv); + if (m_KInvMode) + getKtildeinvFromEcm(E_over_mref,Kv); else - getKtildeinvFromEcm(E_over_mref,Kv);} + getKtildeFromEcm(E_over_mref,Kv);} } @@ -944,7 +841,7 @@ void BoxQuantization::get_ktilde_matrix(double E_over_mref, RealSymmetricMatrix& else{ double kres=evalptr->calculate(rt->getJtimestwo(), rt->getL(),rt->getStimestwo(),rt->getChannelIndex(), - ct->getL(),ct->getStimestwo(),ct->getChannelIndex(),m_kappa_params,Ecm); + ct->getL(),ct->getStimestwo(),ct->getChannelIndex(),Ecm); assign(kres,row,col,herm,Kh,K);} }} } @@ -978,10 +875,10 @@ void BoxQuantization::output_matrices(double E_over_mref, bool Elab, std::ostrea string Kname; if (Elab){ fout << "Elab:="<::iterator rt=m_basis.begin();rt!=m_basis.end();rt++,row++){ uint ap=rt->getChannelIndex(); @@ -1034,7 +931,7 @@ double BoxQuantization::get_determinant(uint N, const RealSymmetricMatrix& Kv, const ComplexHermitianMatrix& B, uint Ndet) { RealDeterminantRoot DR; - if (m_Kinv!=0){ // det( Kinv - B ) + if (m_KInvMode){ // det( Kinv - B ) ComplexHermitianMatrix Q(N); for (uint row=0;row BoxQuantization::get_eigenvalues(double E_over_mref, bool El assign_matrices(E_over_mref,Elab,B,Kv); uint N=B.size(); ComplexHermitianMatrix Q(N); - if (m_Kinv!=0){ // Q = Kinv - B + if (m_KInvMode){ // Q = Kinv - B for (uint row=0;row > m_boxes; std::list m_wzetas; - std::vector m_kappa_params; - KtildeMatrixCalculator* m_Kmat; - KtildeInverseCalculator* m_Kinv; + KtildeMatrixBase* m_Kptr; + bool m_KInvMode; // Prevent copying and no default. @@ -413,26 +412,15 @@ class BoxQuantization public: BoxQuantization(XMLHandler& xmlin, - KtildeMatrixCalculator *Kmatptr); - - BoxQuantization(XMLHandler& xmlin, - KtildeInverseCalculator *Kinvptr); - - BoxQuantization(XMLHandler& xmlin, - KtildeMatrixCalculator *Kmatptr, - KtildeInverseCalculator *Kinvptr); // one of the pointers must be null - - BoxQuantization(const std::string& mom_ray, uint mom_int_sq, - const std::string& lgirrep, - const std::vector& chan_infos, - const std::vector lmaxes, - KtildeMatrixCalculator *Kmatptr); + KtildeMatrixBase *Kptr, + bool KInvMode); BoxQuantization(const std::string& mom_ray, uint mom_int_sq, const std::string& lgirrep, const std::vector& chan_infos, const std::vector lmaxes, - KtildeInverseCalculator *Kinvptr); + KtildeMatrixBase *Kptr, + bool KInvMode); ~BoxQuantization(); @@ -459,8 +447,6 @@ class BoxQuantization void setMassesOverRef(uint channel_index, double mass1_over_ref, double mass2_over_ref); - void setKtildeParameters(std::vector kappa_params); - double getRefMassL() const {return m_mref_L;} @@ -472,14 +458,10 @@ class BoxQuantization {return m_masses2.at(channel_index);} bool isKtildeInverseMode() const - {return (m_Kinv!=0);} + {return m_KInvMode;} bool isKtildeMode() const - {return (m_Kmat!=0);} - - uint getNumberOfKtildeParameters() const - {if (m_Kmat!=0) return m_Kmat->getNumberOfParameters(); - else return m_Kinv->getNumberOfParameters();} + {return !m_KInvMode;} std::string output(int indent=0) const; // XML output @@ -492,19 +474,6 @@ class BoxQuantization std::string outputBasis(int indent=0) const; // XML output - void outputKFitParams(XMLHandler& xmlout) const; // XML output - - std::string outputKFitParams(int indent=0) const; // XML output - - - const std::vector& getFitParameterInfos() const; - - int getParameterIndex(const KFitParamInfo& kinfo) const; // returns -1 if not found - - double getParameterValue(const KFitParamInfo& kinfo) const; - - std::set getElementInfos() const; - double getEcmOverMrefFromElab(double Elab_over_mref) const; diff --git a/source/driver1.cc b/source/driver1.cc index 4d25268..02edc30 100644 --- a/source/driver1.cc +++ b/source/driver1.cc @@ -19,24 +19,35 @@ string fname(argv[1]); XMLHandler xmlin; xmlin.set_from_file(fname); -KtildeMatrixCalculator *Kmat=0; -KtildeInverseCalculator *Kinv=0; +KtildeMatrixBase *Kptr=0; +bool kInvMode; int k1=xmlin.count_among_children("KtildeMatrix"); int k2=xmlin.count_among_children("KtildeMatrixInverse"); if ((k1+k2)!=1) throw(std::invalid_argument("A single KtildeMatrix or KtildeMatrixInverse tag must be present")); if (k1==1){ - Kmat=new KtildeMatrixCalculator(xmlin);} + Kptr=new KtildeMatrixCalculator(xmlin); + kInvMode = false;} else{ - Kinv=new KtildeInverseCalculator(xmlin);} - -if (Kmat!=0) cout <output()<output()< Kpars; +if (kInvMode) { + Kpars.resize(static_cast(Kptr)->getNumberOfParameters()); + static_cast(Kptr)->setKtildeParameters(Kpars); + cout <(Kptr)->output()<(Kptr)->getNumberOfParameters()); + static_cast(Kptr)->setKtildeParameters(Kpars); + cout <(Kptr)->output()< Kpars; +Kpars.resize(Kinv.getNumberOfParameters()); +Kinv.setKtildeParameters(Kpars); + string mom_ray("ar"); uint mom_int_sq=0; string lgirrep("T1u"); @@ -26,7 +30,7 @@ vector lmaxes; chan_infos.push_back(DecayChannelInfo("pion","pion",0,0,false,true)); lmaxes.push_back(5); -BoxQuantization BQ(mom_ray,mom_int_sq,lgirrep,chan_infos,lmaxes,&Kinv); +BoxQuantization BQ(mom_ray,mom_int_sq,lgirrep,chan_infos,lmaxes,&Kinv, true); cout << BQ.output()< Date: Wed, 18 Apr 2018 16:25:07 +0200 Subject: [PATCH 3/4] Change normalization of box matrix elements. This change in normalization induces a change in the normalization of the Ktilde matrix to make it more natural. The difference is most easily seen by comparing eq. (2.9) of https://arxiv.org/abs/1802.03100v1 (new convention) eq. (46) of https://arxiv.org/abs/1707.05817v1 (old convention) In the old normalization, the K matrix depended on m_ref L, i.e. a finite-volume quantity. Using the new normalization, the K matrix does not pick up a dependence on any finite-volume quantity, and the Breit-Wigner parametrization takes its canonical form. --- source/box_quant.cc | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/source/box_quant.cc b/source/box_quant.cc index e70fa59..5d83665 100644 --- a/source/box_quant.cc +++ b/source/box_quant.cc @@ -794,6 +794,8 @@ void BoxQuantization::assign_matrices(double E_over_mref, bool Elab, ComplexHerm void BoxQuantization::get_box_matrix(double E_over_mref, ComplexHermitianMatrix& Bh, CMatrix& B, bool Elab, bool herm) { + const double pi=3.14159265358979323846264; + if (Elab) for (list >::const_iterator it=m_boxes.begin();it!=m_boxes.end();it++){ it->first->setElementsFromElab(E_over_mref);} @@ -815,7 +817,7 @@ void BoxQuantization::get_box_matrix(double E_over_mref, ComplexHermitianMatrix& else{ BoxMatrixQuantumNumbers bqn(rt->getJtimestwo(),rt->getL(),rt->getOccurrence(), ct->getJtimestwo(),ct->getL(),ct->getOccurrence()); - assign(bptr->getElement(bqn),row,col,herm,Bh,B);} + assign((bptr->getElement(bqn))/pow(m_mref_L/(2.0*pi),(rt->getL())+(ct->getL())+1),row,col,herm,Bh,B);} }} } From e20515caf257ca1e97601e03dda68cbeb5a5679d Mon Sep 17 00:00:00 2001 From: Ben Horz Date: Wed, 18 Apr 2018 16:41:09 +0200 Subject: [PATCH 4/4] Add health warning to README about normalization. --- README | 2 ++ 1 file changed, 2 insertions(+) diff --git a/README b/README index 4b04256..e7ecd22 100644 --- a/README +++ b/README @@ -23,6 +23,8 @@ An overview of this software is given in the paper below: SLACcitation = "%%CITATION = ARXIV:1707.05817;%%" } +**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.