diff --git a/Hadrons/Modules/MDistil/BContraction.hpp b/Hadrons/Modules/MDistil/BContraction.hpp index 3e8c523f..62300c3a 100644 --- a/Hadrons/Modules/MDistil/BContraction.hpp +++ b/Hadrons/Modules/MDistil/BContraction.hpp @@ -18,18 +18,17 @@ BEGIN_HADRONS_NAMESPACE ******************************************************************************/ BEGIN_MODULE_NAMESPACE(MDistil) - - // general baryon tensor set based on Eigen tensors and Grid-allocated memory // Dimensions: // 0 - ext - external field (momentum, EM field, ...) - // 1 - str - spin-color structure + // 1 - str - dirac structure // 2 - t - timeslice - // 3 - i - left distillation mode index - // 4 - j - middle distillation mode index - // 5 - k - left distillation mode index + // 3 - s - free spin index + // 4 - i - left distillation mode index + // 5 - j - middle distillation mode index + // 6 - k - left distillation mode index // template - // using BaryonTensorSet = Eigen::TensorMap>; + // using BaryonTensorSet = Eigen::TensorMap>; class BContractionPar: Serializable @@ -74,7 +73,7 @@ protected: /*class BFieldIO: Serializable{ public: - using BaryonTensorSet = Eigen::Tensor; + using BaryonTensorSet = Eigen::Tensor; GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO, BaryonTensorSet, BField ); @@ -183,9 +182,9 @@ void TBContraction::execute(void) Gamma::Algebra::GammaYGamma5, // i gamma_4 C gamma_5 = i gamma_2 gamma_5 }; std::vector factor23{{0.,-1.},{0.,1.},{0.,1.}}; - using BaryonTensorSet = Eigen::Tensor; + using BaryonTensorSet = Eigen::Tensor; int Ngamma=3; - BaryonTensorSet BField3(Nmom,4*Ngamma,Nt,N_1,N_2,N_3); + BaryonTensorSet BField3(Nmom,Ngamma,Nt,4,N_1,N_2,N_3); Eigen::Tensor corr(Nmom,4,Nt); @@ -216,7 +215,7 @@ void TBContraction::execute(void) tmp111 = 0.5*(double)parity*(tmp111 + tmp222); // P_\pm * ... diquark2 = factor23[0]*innerProduct(tmp22s,tmp333); for (int is=0 ; is < 4 ; is++){ - BField3(imom,is+4*ig,t,i1,i2,i3)+=(double)epsilon_sgn[ie]*tmp111()(is)()*diquark2; + BField3(imom,ig,t,is,i1,i2,i3)+=(double)epsilon_sgn[ie]*tmp111()(is)()*diquark2; } } } @@ -228,7 +227,7 @@ void TBContraction::execute(void) } for (int is=0 ; is < 4 ; is++){ for (int t=0 ; t < Nt ; t++){ - std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField3(0,is,t,0,0,0) << std::endl; + std::cout << "BaryonField(is=" << is << ",t=" << t << ") = " << BField3(0,0,t,is,0,0,0) << std::endl; } } @@ -240,45 +239,6 @@ void TBContraction::execute(void) Hdf5Writer writer(filename); write(writer,"BaryonField",BField_save.BField); - /* int Npairs = 0; - char left[] = "uud"; - char right[] = "uud"; - std::vector pairs(6); - for (int ie=0, i=0 ; ie < 6 ; ie++){ - if (left[0] == right[epsilon[ie][0]] && left[1] == right[epsilon[ie][1]] && left[2] == right[epsilon[ie][2]]){ - pairs[i] = ie; - i++; - Npairs++; - } - } - pairs.resize(Npairs); - std::cout << Npairs << " pairs: " << pairs << std::endl; - for (int imom=0 ; imom < Nmom ; imom++){ - for (int is=0 ; is < 4 ; is++){ - for (int t=0 ; t < Nt ; t++){ - corr(imom,is,t) = 0.; - } - } - } - for (int ipair=0 ; ipair < Npairs ; ipair++){ - Eigen::array, 3> product_dims = { Eigen::IndexPair(0,epsilon[pairs[ipair]][0]),Eigen::IndexPair(1,epsilon[pairs[ipair]][1]) ,Eigen::IndexPair(2,epsilon[pairs[ipair]][2]) }; - for (int imom=0 ; imom < Nmom ; imom++){ - Eigen::Tensor B5 = BField3.chip(imom,0); - for (int is=0 ; is < 4 ; is++){ - Eigen::Tensor B4 = B5.chip(is,0); - for (int t=0 ; t < Nt ; t++){ - Eigen::Tensor B3 = B4.chip(t,0); - Eigen::Tensor C2 = B3.contract(B3,product_dims); - corr(imom,is,t) += (double)epsilon_sgn[pairs[ipair]]*C2(0); - } - } - } - } - for (int is=0 ; is < 4 ; is++){ - for (int t=0 ; t < Nt ; t++){ - std::cout << "C2(is=" << is << ",t=" << t << ") = " << corr(0,is,t) << std::endl; - } - }*/ } END_MODULE_NAMESPACE diff --git a/Hadrons/Modules/MDistil/Baryon2pt.hpp b/Hadrons/Modules/MDistil/Baryon2pt.hpp index e97dab16..75a23e7a 100644 --- a/Hadrons/Modules/MDistil/Baryon2pt.hpp +++ b/Hadrons/Modules/MDistil/Baryon2pt.hpp @@ -24,6 +24,8 @@ public: GRID_SERIALIZABLE_CLASS_MEMBERS(Baryon2ptPar, std::string, inputL, std::string, inputR, + std::string, quarksL, + std::string, quarksR, std::string, output ); }; @@ -95,6 +97,8 @@ void TBaryon2pt::execute(void) const std::string &inputL{par().inputL}; const std::string &inputR{par().inputR}; + const std::string &inputL{par().quarksL}; + const std::string &inputR{par().quarksR}; const std::string &output{par().output}; int Nmom=1; @@ -108,10 +112,10 @@ void TBaryon2pt::execute(void) int N_2=20; int N_3=20; - // using BaryonTensorSet = Eigen::Tensor; + // using BaryonTensorSet = Eigen::Tensor; BFieldIO BFieldL; - BFieldL.BField.resize(Nmom,4*Ngamma,Nt,N_1,N_2,N_3); + BFieldL.BField.resize(Nmom,Ngamma,Nt,4,N_1,N_2,N_3); std::string filenameL ="./" + inputL + ".h5"; std::cout << "Reading from file " << filenameL << std::endl; @@ -119,7 +123,7 @@ void TBaryon2pt::execute(void) read(readerL,"BaryonField",BFieldL.BField); BFieldIO BFieldR; - BFieldR.BField.resize(Nmom,4*Ngamma,Nt,N_1,N_2,N_3); + BFieldR.BField.resize(Nmom,Ngamma,Nt,4,N_1,N_2,N_3); std::string filenameR ="./" + inputR + ".h5"; std::cout << "Reading from file " << filenameR << std::endl; @@ -153,16 +157,20 @@ void TBaryon2pt::execute(void) Eigen::array, 3> product_dims = { Eigen::IndexPair(0,epsilon[pairs[ipair]][0]),Eigen::IndexPair(1,epsilon[pairs[ipair]][1]) ,Eigen::IndexPair(2,epsilon[pairs[ipair]][2]) }; for (int imom=0 ; imom < Nmom ; imom++){ std::cout << imom << std::endl; - Eigen::Tensor B5L = BFieldL.BField.chip(imom,0); - Eigen::Tensor B5R = BFieldR.BField.chip(imom,0); - for (int is=0 ; is < 4 ; is++){ - Eigen::Tensor B4L = B5L.chip(is,0); - Eigen::Tensor B4R = B5R.chip(is,0); + Eigen::Tensor B6L = BFieldL.BField.chip(imom,0); + Eigen::Tensor B6R = BFieldR.BField.chip(imom,0); + for (int ig=0 ; ig < Ngamma ; ig++){ + Eigen::Tensor B5L = B6L.chip(ig,0); + Eigen::Tensor B5R = B6R.chip(ig,0); for (int t=0 ; t < Nt ; t++){ - Eigen::Tensor B3L = B4L.chip(t,0); - Eigen::Tensor B3R = B4R.chip(tsrc,0); - Eigen::Tensor C2 = B3L.contract(B3R,product_dims); - corr(imom,t) += (double)epsilon_sgn[pairs[ipair]]*C2(0); + Eigen::Tensor B4L = B5L.chip(t,0); + Eigen::Tensor B4R = B5R.chip(tsrc,0); + for (int is=0 ; is < 4 ; is++){ + Eigen::Tensor B3L = B4L.chip(is,0); + Eigen::Tensor B3R = B4R.chip(is,0); + Eigen::Tensor C2 = B3L.contract(B3R,product_dims); + corr(imom,t) += (double)epsilon_sgn[pairs[ipair]]*C2(0); + } } } } diff --git a/Hadrons/Modules/MDistil/Distil.hpp b/Hadrons/Modules/MDistil/Distil.hpp index d6673dd0..af7e2786 100644 --- a/Hadrons/Modules/MDistil/Distil.hpp +++ b/Hadrons/Modules/MDistil/Distil.hpp @@ -202,7 +202,7 @@ public: class BFieldIO: Serializable{ public: - using BaryonTensorSet = Eigen::Tensor; + using BaryonTensorSet = Eigen::Tensor; GRID_SERIALIZABLE_CLASS_MEMBERS(BFieldIO, BaryonTensorSet, BField ); diff --git a/tests/hadrons/Test_hadrons_distil.cc b/tests/hadrons/Test_hadrons_distil.cc index bfa8c305..5b43fac0 100644 --- a/tests/hadrons/Test_hadrons_distil.cc +++ b/tests/hadrons/Test_hadrons_distil.cc @@ -292,6 +292,8 @@ void test_Baryon2pt(Application &application) MDistil::Baryon2pt::Par Baryon2ptPar; Baryon2ptPar.inputL="BaryonFieldPhi"; Baryon2ptPar.inputR="BaryonFieldRho"; + Baryon2ptPar.quarksL="uud"; + Baryon2ptPar.quarksR="uud"; Baryon2ptPar.output="C2_baryon"; application.createModule("C2_b",Baryon2ptPar); }