From 868c33e6512844908fe6c906f33cf0af3be48525 Mon Sep 17 00:00:00 2001 From: Alexander Schmidt Date: Tue, 23 Dec 2014 11:29:24 +0100 Subject: [PATCH] Added missing isospin correlators. --- u1casc-ordinary/obs_corriso00.hpp | 150 ++++++++++++++++++++++++++++ u1casc-ordinary/obs_corriso10.hpp | 150 ++++++++++++++++++++++++++++ u1casc-ordinary/u1casc-ordinary.cpp | 10 ++ 3 files changed, 310 insertions(+) create mode 100644 u1casc-ordinary/obs_corriso00.hpp create mode 100644 u1casc-ordinary/obs_corriso10.hpp diff --git a/u1casc-ordinary/obs_corriso00.hpp b/u1casc-ordinary/obs_corriso00.hpp new file mode 100644 index 0000000..ca4f2db --- /dev/null +++ b/u1casc-ordinary/obs_corriso00.hpp @@ -0,0 +1,150 @@ +#ifndef OBS_CORRISO00_HPP +#define OBS_CORRISO00_HPP + +#include "latlib/o815/o815.h" + +#include "latlib/writeout.h" + +#include "latlib/obstat.hpp" + +#include +#include + +#include + +using namespace std; + +class obs_corriso00 : public o815::obs { + +public: + obs_corriso00(o815 *_O815); + +private: + void _start(); + void _meas(bool loadedobs, const int& nthmeas); + void _finish(); + + void corrCompute(); + static complex effMass(vector < complex > *preCalculated, vector< complex > *excludedMeas, int nmeas, void *para); + static void preEffMass(vector< vector < complex > > *allVals, vector < complex > *preCalculated, void *para); + + sim *Sim; + obstat< complex, complex > oC; + + int spatialV; + + complex *OM; +}; + +obs_corriso00::obs_corriso00(o815 *_O815) : o815::obs("corriso00", + _O815->paraQ->getParaNames() + + "tsep" + ":iso00_real:iso00_imag:iso00_abs:iso00_mass:iso00_mass_err", + _O815, sizeof(complex)*(_O815->comargs.lsize[1]/2+1) ) { + + OM = (complex*)(obsMem); + + Sim = (sim*)O815->Sim; + spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0]; +} + +void obs_corriso00::_start() { + //*out << "OBS_test: start" << endl; +}; + +void obs_corriso00::_meas(bool loadedobs, const int& nthmeas) { + if (!loadedobs) + corrCompute(); + + oC.addMeas( OM, O815->comargs.lsize[1]/2+1 ); +}; + +void obs_corriso00::_finish() { + int compid_corr[O815->comargs.lsize[1]/2]; + int compid_effmass[O815->comargs.lsize[1]/2-1]; + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) + compid_corr[itsep] = oC.computeEasy(itsep); + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) { + pair effmasspass( itsep, O815->comargs.lsize[1]/2 ); + compid_effmass[itsep] = oC.computeJack(obs_corriso00::preEffMass, obs_corriso00::effMass, &effmasspass); + } + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) { + *out << O815->paraQ->getParaVals(); + *out << "\t" << itsep; + + *out << "\t" << real( oC.getMean(compid_corr[itsep]) ) + << "\t" << imag( oC.getMean(compid_corr[itsep]) ) + << "\t" << abs ( oC.getMean(compid_corr[itsep]) ); + + if ( itsep < O815->comargs.lsize[1]/2-1 ) + *out << "\t" << real( oC.getMean(compid_effmass[itsep]) ) + << "\t" << real( oC.getErr (compid_effmass[itsep]) ); + else + *out << "\t" << NAN << "\t" << NAN; + + *out << endl; + } + + oC.reset(); +} + + +void obs_corriso00::corrCompute() +{ + complex phislice[O815->comargs.lsize[1]]; + + OM[O815->comargs.lsize[1]/2] = 0; + + for (int it = 0; it < O815->comargs.lsize[1]; it++) { + phislice[it] = 0; + + for (int ix = 0; ix < spatialV; ix++) + phislice[it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ] + + conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]; + + phislice[it] /= sqrt(2)*spatialV; + + OM[O815->comargs.lsize[1]/2] += phislice[it]; + } + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) { + OM[itsep] = 0; + + for (int it = 0; it < O815->comargs.lsize[1]; it++) + OM[itsep] += phislice[ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[it] ); + + OM[itsep] /= O815->comargs.lsize[1]; + } + + OM[O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1]; +} + +void obs_corriso00::preEffMass(vector< vector < complex > > *allVals, vector < complex > *preCalculated, void *para) { + pair *myparas = (pair*)para; + + preCalculated->push_back(0); + preCalculated->push_back(0); + preCalculated->push_back(0); + + for(vector< vector< complex > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) { + (*preCalculated)[0] += (*valIt)[myparas->first]; + (*preCalculated)[1] += (*valIt)[myparas->first+1]; + (*preCalculated)[2] += (*valIt)[myparas->second]; + } +} + +complex obs_corriso00::effMass(vector < complex > *preCalculated, vector< complex > *excludedMeas, int nmeas, void *para) { + pair *myparas = (pair*)para; + + double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex)(nmeas-1) ); + + return std::log( abs( + ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex)(nmeas-1) - disconnected ) / + ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex)(nmeas-1) - disconnected ) + ) ); +} + +#endif diff --git a/u1casc-ordinary/obs_corriso10.hpp b/u1casc-ordinary/obs_corriso10.hpp new file mode 100644 index 0000000..3bdb7b5 --- /dev/null +++ b/u1casc-ordinary/obs_corriso10.hpp @@ -0,0 +1,150 @@ +#ifndef OBS_CORRISO10_HPP +#define OBS_CORRISO10_HPP + +#include "latlib/o815/o815.h" + +#include "latlib/writeout.h" + +#include "latlib/obstat.hpp" + +#include +#include + +#include + +using namespace std; + +class obs_corriso10 : public o815::obs { + +public: + obs_corriso10(o815 *_O815); + +private: + void _start(); + void _meas(bool loadedobs, const int& nthmeas); + void _finish(); + + void corrCompute(); + static complex effMass(vector < complex > *preCalculated, vector< complex > *excludedMeas, int nmeas, void *para); + static void preEffMass(vector< vector < complex > > *allVals, vector < complex > *preCalculated, void *para); + + sim *Sim; + obstat< complex, complex > oC; + + int spatialV; + + complex *OM; +}; + +obs_corriso10::obs_corriso10(o815 *_O815) : o815::obs("corriso10", + _O815->paraQ->getParaNames() + + "tsep" + ":iso10_real:iso10_imag:iso10_abs:iso10_mass:iso10_mass_err", + _O815, sizeof(complex)*(_O815->comargs.lsize[1]/2+1) ) { + + OM = (complex*)(obsMem); + + Sim = (sim*)O815->Sim; + spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0]; +} + +void obs_corriso10::_start() { + //*out << "OBS_test: start" << endl; +}; + +void obs_corriso10::_meas(bool loadedobs, const int& nthmeas) { + if (!loadedobs) + corrCompute(); + + oC.addMeas( OM, O815->comargs.lsize[1]/2+1 ); +}; + +void obs_corriso10::_finish() { + int compid_corr[O815->comargs.lsize[1]/2]; + int compid_effmass[O815->comargs.lsize[1]/2-1]; + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) + compid_corr[itsep] = oC.computeEasy(itsep); + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) { + pair effmasspass( itsep, O815->comargs.lsize[1]/2 ); + compid_effmass[itsep] = oC.computeJack(obs_corriso10::preEffMass, obs_corriso10::effMass, &effmasspass); + } + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) { + *out << O815->paraQ->getParaVals(); + *out << "\t" << itsep; + + *out << "\t" << real( oC.getMean(compid_corr[itsep]) ) + << "\t" << imag( oC.getMean(compid_corr[itsep]) ) + << "\t" << abs ( oC.getMean(compid_corr[itsep]) ); + + if ( itsep < O815->comargs.lsize[1]/2-1 ) + *out << "\t" << real( oC.getMean(compid_effmass[itsep]) ) + << "\t" << real( oC.getErr (compid_effmass[itsep]) ); + else + *out << "\t" << NAN << "\t" << NAN; + + *out << endl; + } + + oC.reset(); +} + + +void obs_corriso10::corrCompute() +{ + complex phislice[O815->comargs.lsize[1]]; + + OM[O815->comargs.lsize[1]/2] = 0; + + for (int it = 0; it < O815->comargs.lsize[1]; it++) { + phislice[it] = 0; + + for (int ix = 0; ix < spatialV; ix++) + phislice[it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ] + - conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]; + + phislice[it] /= sqrt(2)*spatialV; + + OM[O815->comargs.lsize[1]/2] += phislice[it]; + } + + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) { + OM[itsep] = 0; + + for (int it = 0; it < O815->comargs.lsize[1]; it++) + OM[itsep] += phislice[ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[it] ); + + OM[itsep] /= O815->comargs.lsize[1]; + } + + OM[O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1]; +} + +void obs_corriso10::preEffMass(vector< vector < complex > > *allVals, vector < complex > *preCalculated, void *para) { + pair *myparas = (pair*)para; + + preCalculated->push_back(0); + preCalculated->push_back(0); + preCalculated->push_back(0); + + for(vector< vector< complex > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) { + (*preCalculated)[0] += (*valIt)[myparas->first]; + (*preCalculated)[1] += (*valIt)[myparas->first+1]; + (*preCalculated)[2] += (*valIt)[myparas->second]; + } +} + +complex obs_corriso10::effMass(vector < complex > *preCalculated, vector< complex > *excludedMeas, int nmeas, void *para) { + pair *myparas = (pair*)para; + + double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex)(nmeas-1) ); + + return std::log( abs( + ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex)(nmeas-1) - disconnected ) / + ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex)(nmeas-1) - disconnected ) + ) ); +} + +#endif diff --git a/u1casc-ordinary/u1casc-ordinary.cpp b/u1casc-ordinary/u1casc-ordinary.cpp index 4803fe1..b6b69c0 100644 --- a/u1casc-ordinary/u1casc-ordinary.cpp +++ b/u1casc-ordinary/u1casc-ordinary.cpp @@ -17,6 +17,8 @@ sim *Sim; #include "obs_corrphichi.hpp" #include "obs_corrphi.hpp" #include "obs_corrchi.hpp" +#include "obs_corriso10.hpp" +#include "obs_corriso00.hpp" o815::comoption specOps[] = { { "kappaone", required_argument, NULL, 'r', "set inverse mass kappa_1", "min:max:inc" }, @@ -95,6 +97,14 @@ void parseLonelyArgs() *O815->out->log << "MASTER: registered observable: corrchi" << endl << flush; O815->observables.push_back(new obs_corrchi(O815)); } + else if ( strcmp(*lonit, "corriso10") == 0 ) { + *O815->out->log << "MASTER: registered observable: corriso10" << endl << flush; + O815->observables.push_back(new obs_corriso10(O815)); + } + else if ( strcmp(*lonit, "corriso00") == 0 ) { + *O815->out->log << "MASTER: registered observable: corriso00" << endl << flush; + O815->observables.push_back(new obs_corriso00(O815)); + } } } -- 2.39.5