From 83f20c4a1509e09d4836e2e15fc8e142950e15ec Mon Sep 17 00:00:00 2001 From: Alexander Schmidt Date: Thu, 27 Feb 2014 23:06:23 +0100 Subject: [PATCH] Migrated observable plaq. --- u1casc-flux/obs_plaq.hpp | 126 ++++++++++++++++++++++++++++++++++++ u1casc-flux/sim.hpp | 4 +- u1casc-flux/u1casc-flux.cpp | 10 +-- 3 files changed, 134 insertions(+), 6 deletions(-) create mode 100644 u1casc-flux/obs_plaq.hpp diff --git a/u1casc-flux/obs_plaq.hpp b/u1casc-flux/obs_plaq.hpp new file mode 100644 index 0000000..249f091 --- /dev/null +++ b/u1casc-flux/obs_plaq.hpp @@ -0,0 +1,126 @@ +#ifndef OBS_PLAQ_HPP +#define OBS_PLAQ_HPP + +#include "latlib/o815/o815.h" + +#include "latlib/writeout.h" + +#include "latlib/obstat.hpp" + +#include + +using namespace std; + +class obs_plaq : public o815::obs { + +public: + struct obsmem { + double A; + double B; + }; + obs_plaq(o815 *_O815); + +private: + void _start(); + void _meas(bool loadedobs, const int& nthmeas); + void _finish(); + + obsmem* OM; + + static void prePlaqSus(vector< vector > *allVals, vector *preCalculated, void *para); + static double plaqSus(vector *preCalculated, vector *excludedMeas, int nmeas, void *para); + + sim *Sim; + obstat oPlaq; + + double A(const double& beta); + double B(const double& beta); +}; + +obs_plaq::obs_plaq(o815 *_O815) : o815::obs("plaq", + _O815->paraQ->getParaNames() + "plaq:plaq_err:plaqsus:plaqsus_err", + _O815, sizeof(obsmem) ) { + OM = (obsmem*)obsMem; + Sim = (sim*)O815->Sim; +} + +void obs_plaq::_start() { + //*out << "OBS_test: start" << endl; +}; + +void obs_plaq::_meas(bool loadedobs, const int& nthmeas) { + if (!loadedobs) { + OM->A = A((*paraQ)["beta"]); + OM->B = B((*paraQ)["beta"]); + } + + double tmpMeas[2] = {OM->A, OM->B}; + oPlaq.addMeas(tmpMeas, 2); +}; + +void obs_plaq::_finish() { + *out << O815->paraQ->getParaVals(); + + int compid_plaq, compid_plaqsus; + + compid_plaq = oPlaq.computeEasy(); + compid_plaqsus = oPlaq.computeJack(obs_plaq::prePlaqSus, obs_plaq::plaqSus, &(Sim->lsize4)); + + *out << "\t" << oPlaq.getMean(compid_plaq) << "\t" << oPlaq.getErr(compid_plaq); + *out << "\t" << oPlaq.getMean(compid_plaqsus) << "\t" << oPlaq.getErr(compid_plaqsus); + + *out << endl; + + oPlaq.reset(); +}; + +void obs_plaq::prePlaqSus(vector< vector > *allVals, vector *preCalculated, void *para) { + preCalculated->push_back(0); + preCalculated->push_back(0); + + for(vector< vector >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) { + (*preCalculated)[0] += (*valIt)[1] + pow((*valIt)[0], 2)*(*(int*)para*6); + (*preCalculated)[1] += (*valIt)[0]; + } +} + +double obs_plaq::plaqSus(vector *preCalculated, vector *excludedMeas, int nmeas, void *para) { + return ( (*preCalculated)[0] - pow((*excludedMeas)[0], 2) * *(int*)para * 6 - (*excludedMeas)[1] ) / (nmeas-1) + - pow( ( (*preCalculated)[1] - (*excludedMeas)[0] ) / (nmeas-1), 2 ) + * *(int*)para * 6; +} + +double obs_plaq::A(const double& beta) +{ + double A = 0; + + for( int x=0; xlsize4; x++ ) + for( int sigma=0; sigma<3; sigma++ ) + for( int tau=sigma+1; tau<4; tau++ ) { + A += Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] ); + A += Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] / beta; + } + + return A/(6*Sim->lsize4); +} + +double obs_plaq::B(const double& beta) +{ + double B=0; + + for( int x=0; xlsize4; x++ ) + for( int sigma=0; sigma<3; sigma++ ) + for( int tau=sigma+1; tau<4; tau++ ) { + B += Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 2 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] ); + B += ( (Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) ) / + ( beta * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ]) ); + B -= pow(Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1 ) / Sim->I( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] ), 2); + B -= ( Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ] + 1) ) / + ( beta * Sim->I(Sim->p[ x*6 + Sim->iPlaq(sigma,tau) ]) ); + B -= Sim->p[ x*6 + Sim->iPlaq(sigma,tau)] / pow(beta,2); + } + + return B/(6*Sim->lsize4); +} + +#endif diff --git a/u1casc-flux/sim.hpp b/u1casc-flux/sim.hpp index 0da2b8a..7fd59af 100644 --- a/u1casc-flux/sim.hpp +++ b/u1casc-flux/sim.hpp @@ -22,6 +22,8 @@ public: double kappa[2], lambda[2], Mu[2], beta; int wArg(int x, int i); double WF(int nom, int denom, int iflav); + int iPlaq(int i, int j); + double I(int p); private: void _makeSweep(); @@ -32,9 +34,7 @@ private: int cubeUpdate(int x, int orient); void getCube(int cube[3], int orient); int plaqUpdate(int x, int mu, int nu, int i); - int iPlaq(int i, int j); int lpUpdate(int x, int mu, int i); - double I(int p); static double wkern(double t, void *params); double wGsl(int n, int iflav); struct bcache { diff --git a/u1casc-flux/u1casc-flux.cpp b/u1casc-flux/u1casc-flux.cpp index d0d6347..64869e4 100644 --- a/u1casc-flux/u1casc-flux.cpp +++ b/u1casc-flux/u1casc-flux.cpp @@ -9,7 +9,7 @@ o815 *O815; sim *Sim; #include "obs_phi2.hpp" -//#include "obs_plaq.hpp" +#include "obs_plaq.hpp" //#include "obs_corr.hpp" o815::comoption specOps[] = { @@ -77,15 +77,17 @@ void parseLonelyArgs() if ( strcmp(*lonit, "phi2") == 0 ) { *O815->out->log << "MASTER: registered observable: phi2" << endl << flush; O815->observables.push_back(new obs_phi2(O815)); - } /* + } else if ( strcmp(*lonit, "plaq") == 0 ) { *O815->out->log << "MASTER: registered observable: plaq" << endl << flush; O815->observables.push_back(new obs_plaq(O815)); } - else if ( strcmp(*lonit, "corr") == 0 ) { + /* + else if ( strcmp(*lonit, "corr") == 0 ) { *O815->out->log << "MASTER: registered observable: corr" << endl << flush; O815->observables.push_back(new obs_corr(O815)); - } */ + } + */ } } -- 2.39.5