+++ /dev/null
-#ifndef OBS_CORR_HPP
-#define OBS_CORR_HPP
-
-#include "latlib/o815/o815.h"
-
-#include "latlib/writeout.h"
-
-#include "latlib/obstat.hpp"
-
-#include <iostream>
-#include <complex>
-
-#include <cmath>
-
-#define U1 0
-#define U2 1
-#define F1 2
-#define F2 3
-
-using namespace std;
-
-class obs_corr : public o815::obs {
-
-public:
- obs_corr(o815 *_O815);
-
-private:
- void _start();
- void _meas(bool loadedobs, const int& nthmeas);
- void _finish();
-
- void corrCompute();
- static complex<double> effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para);
- static void preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para);
-
- sim *Sim;
- obstat< complex<double>, complex<double> > oC[4];
- obstat< complex<double>, complex<double> > oD[4];
-
- int spatialV;
-
- complex<double> *OM[4];
-};
-
-obs_corr::obs_corr(o815 *_O815) : o815::obs("corr",
- _O815->paraQ->getParaNames() +
- "tsep"
- ":u1_real:u1_imag:u1_abs:u1_effmass:u1_effmass_err"
- ":u2_real:u2_imag:u2_abs:u2_effmass:u2_effmass_err"
- ":f1_real:f1_imag:f1_abs:f1_effmass:f1_effmass_err"
- ":f2_real:f2_imag:f2_abs:f2_effmass:f2_effmass_err",
- _O815, 4*sizeof(complex<double>)*(_O815->comargs.lsize[1]/2+1) ) {
- for (int ivar = 0; ivar<4; ivar++)
- OM[ivar] = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 + 1 ) );
-
- Sim = (sim*)O815->Sim;
- spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
-}
-
-void obs_corr::_start() {
- //*out << "OBS_test: start" << endl;
-};
-
-void obs_corr::_meas(bool loadedobs, const int& nthmeas) {
- if (!loadedobs)
- corrCompute();
-
- oC[U1].addMeas( OM[U1], O815->comargs.lsize[1]/2+1 );
- oC[U2].addMeas( OM[U2], O815->comargs.lsize[1]/2+1 );
- oC[F1].addMeas( OM[F1], O815->comargs.lsize[1]/2+1 );
- oC[F2].addMeas( OM[F2], O815->comargs.lsize[1]/2+1 );
-};
-
-void obs_corr::_finish() {
- int compid_corr[O815->comargs.lsize[1]/2][4];
- int compid_effmass[O815->comargs.lsize[1]/2-1][4];
-
- for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++)
- for (int icorr = 0; icorr < 4; icorr++)
- compid_corr[itsep][icorr] = oC[icorr].computeEasy(itsep);
-
-
- for (int icorr = 0; icorr < 4; icorr++) {
- for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) {
- pair<int,int> effmasspass( itsep, O815->comargs.lsize[1]/2 );
- compid_effmass[itsep][icorr] = oC[icorr].computeJack(obs_corr::preEffMass, obs_corr::effMass, &effmasspass);
- }
- }
-
- for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
- *out << O815->paraQ->getParaVals();
- *out << "\t" << itsep;
- for (int icorr = 0; icorr < 4; icorr++) {
- *out << "\t" << real( oC[icorr].getMean(compid_corr[itsep][icorr]) )
- << "\t" << imag( oC[icorr].getMean(compid_corr[itsep][icorr]) )
- << "\t" << abs ( oC[icorr].getMean(compid_corr[itsep][icorr]) );
-
-
- if ( itsep < O815->comargs.lsize[1]/2-1 )
- *out << "\t" << real( oC[icorr].getMean(compid_effmass[itsep][icorr]) )
- << "\t" << real( oC[icorr].getErr (compid_effmass[itsep][icorr]) );
- else
- *out << "\t" << NAN << "\t" << NAN;
- }
- *out << endl;
- }
-
- for (int icorr = 0; icorr < 4; icorr++) {
- oC[icorr].reset();
- oD[icorr].reset();
- }
-};
-
-void obs_corr::corrCompute()
-{
- complex<double> phislice[4][O815->comargs.lsize[1]];
-
- OM[U1][O815->comargs.lsize[1]/2] = 0;
- OM[U2][O815->comargs.lsize[1]/2] = 0;
- OM[F1][O815->comargs.lsize[1]/2] = 0;
- OM[F2][O815->comargs.lsize[1]/2] = 0;
-
- for (int it = 0; it < O815->comargs.lsize[1]; it++) {
- phislice[U1][it] = 0;
- phislice[U2][it] = 0;
- phislice[F1][it] = 0;
- phislice[F2][it] = 0;
-
- for (int ix = 0; ix < spatialV; ix++) {
- phislice[U1][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
- phislice[U2][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
- phislice[F1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
- phislice[F2][it] += Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
- }
- phislice[U1][it] /= spatialV;
- phislice[U2][it] /= spatialV;
- phislice[F1][it] /= spatialV;
- phislice[F2][it] /= spatialV;
-
- OM[U1][O815->comargs.lsize[1]/2] += phislice[U1][it];
- OM[U2][O815->comargs.lsize[1]/2] += phislice[U2][it];
- OM[F1][O815->comargs.lsize[1]/2] += phislice[F1][it];
- OM[F2][O815->comargs.lsize[1]/2] += phislice[F2][it];
- }
-
- for (int icorr = 0; icorr < 4; icorr++) {
- for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
- OM[icorr][itsep] = 0;
-
- for (int it = 0; it < O815->comargs.lsize[1]; it++)
- OM[icorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[icorr][it] );
-
- OM[icorr][itsep] /= O815->comargs.lsize[1];
- }
-
- OM[icorr][O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1];
- }
-}
-
-void obs_corr::preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para) {
- pair<int,int> *myparas = (pair<int,int>*)para;
-
- preCalculated->push_back(0);
- preCalculated->push_back(0);
- preCalculated->push_back(0);
-
- for(vector< vector< complex<double> > >::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<double> obs_corr::effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para) {
- pair<int,int> *myparas = (pair<int,int>*)para;
-
- double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex<double>)(nmeas-1) );
-
- return std::log( abs(
- ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex<double>)(nmeas-1) - disconnected ) /
- ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex<double>)(nmeas-1) - disconnected )
- ) );
-}
-
-#endif
projects / phys / u1casc.git / blobdiff