4 #include "latlib/o815/o815.h"
6 #include "latlib/writeout.h"
8 #include "latlib/obstat.hpp"
22 class obs_corr : public o815::obs {
25 obs_corr(o815 *_O815);
29 void _meas(bool loadedobs, const int& nthmeas);
33 static complex<double> effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para);
34 static void preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para);
37 obstat< complex<double>, complex<double> > oC[4];
38 obstat< complex<double>, complex<double> > oD[4];
42 complex<double> *OM[4];
45 obs_corr::obs_corr(o815 *_O815) : o815::obs("corr",
46 _O815->paraQ->getParaNames() +
48 ":u1_real:u1_imag:u1_abs:u1_effmass:u1_effmass_err"
49 ":u2_real:u2_imag:u2_abs:u2_effmass:u2_effmass_err"
50 ":f1_real:f1_imag:f1_abs:f1_effmass:f1_effmass_err"
51 ":f2_real:f2_imag:f2_abs:f2_effmass:f2_effmass_err",
52 _O815, 4*sizeof(complex<double>)*(_O815->comargs.lsize[1]/2+1) ) {
53 for (int ivar = 0; ivar<4; ivar++)
54 OM[ivar] = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 + 1 ) );
56 Sim = (sim*)O815->Sim;
57 spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
60 void obs_corr::_start() {
61 //*out << "OBS_test: start" << endl;
64 void obs_corr::_meas(bool loadedobs, const int& nthmeas) {
68 oC[U1].addMeas( OM[U1], O815->comargs.lsize[1]/2+1 );
69 oC[U2].addMeas( OM[U2], O815->comargs.lsize[1]/2+1 );
70 oC[F1].addMeas( OM[F1], O815->comargs.lsize[1]/2+1 );
71 oC[F2].addMeas( OM[F2], O815->comargs.lsize[1]/2+1 );
74 void obs_corr::_finish() {
75 int compid_corr[O815->comargs.lsize[1]/2][4];
76 int compid_effmass[O815->comargs.lsize[1]/2-1][4];
78 for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++)
79 for (int icorr = 0; icorr < 4; icorr++)
80 compid_corr[itsep][icorr] = oC[icorr].computeEasy(itsep);
83 for (int icorr = 0; icorr < 4; icorr++) {
84 for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) {
85 pair<int,int> effmasspass( itsep, O815->comargs.lsize[1]/2 );
86 compid_effmass[itsep][icorr] = oC[icorr].computeJack(obs_corr::preEffMass, obs_corr::effMass, &effmasspass);
90 for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
91 *out << O815->paraQ->getParaVals();
92 *out << "\t" << itsep;
93 for (int icorr = 0; icorr < 4; icorr++) {
94 *out << "\t" << real( oC[icorr].getMean(compid_corr[itsep][icorr]) )
95 << "\t" << imag( oC[icorr].getMean(compid_corr[itsep][icorr]) )
96 << "\t" << abs ( oC[icorr].getMean(compid_corr[itsep][icorr]) );
99 if ( itsep < O815->comargs.lsize[1]/2-1 )
100 *out << "\t" << real( oC[icorr].getMean(compid_effmass[itsep][icorr]) )
101 << "\t" << real( oC[icorr].getErr (compid_effmass[itsep][icorr]) );
103 *out << "\t" << NAN << "\t" << NAN;
108 for (int icorr = 0; icorr < 4; icorr++) {
114 void obs_corr::corrCompute()
116 complex<double> phislice[4][O815->comargs.lsize[1]];
118 OM[U1][O815->comargs.lsize[1]/2] = 0;
119 OM[U2][O815->comargs.lsize[1]/2] = 0;
120 OM[F1][O815->comargs.lsize[1]/2] = 0;
121 OM[F2][O815->comargs.lsize[1]/2] = 0;
123 for (int it = 0; it < O815->comargs.lsize[1]; it++) {
124 phislice[U1][it] = 0;
125 phislice[U2][it] = 0;
126 phislice[F1][it] = 0;
127 phislice[F2][it] = 0;
129 for (int ix = 0; ix < spatialV; ix++) {
130 phislice[U1][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
131 phislice[U2][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
132 phislice[F1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
133 phislice[F2][it] += Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
135 phislice[U1][it] /= spatialV;
136 phislice[U2][it] /= spatialV;
137 phislice[F1][it] /= spatialV;
138 phislice[F2][it] /= spatialV;
140 OM[U1][O815->comargs.lsize[1]/2] += phislice[U1][it];
141 OM[U2][O815->comargs.lsize[1]/2] += phislice[U2][it];
142 OM[F1][O815->comargs.lsize[1]/2] += phislice[F1][it];
143 OM[F2][O815->comargs.lsize[1]/2] += phislice[F2][it];
146 for (int icorr = 0; icorr < 4; icorr++) {
147 for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
148 OM[icorr][itsep] = 0;
150 for (int it = 0; it < O815->comargs.lsize[1]; it++)
151 OM[icorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[icorr][it] );
153 OM[icorr][itsep] /= O815->comargs.lsize[1];
156 OM[icorr][O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1];
160 void obs_corr::preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para) {
161 pair<int,int> *myparas = (pair<int,int>*)para;
163 preCalculated->push_back(0);
164 preCalculated->push_back(0);
165 preCalculated->push_back(0);
167 for(vector< vector< complex<double> > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
168 (*preCalculated)[0] += (*valIt)[myparas->first];
169 (*preCalculated)[1] += (*valIt)[myparas->first+1];
170 (*preCalculated)[2] += (*valIt)[myparas->second];
174 complex<double> obs_corr::effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para) {
175 pair<int,int> *myparas = (pair<int,int>*)para;
177 double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex<double>)(nmeas-1) );
179 return std::log( abs(
180 ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex<double>)(nmeas-1) - disconnected ) /
181 ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex<double>)(nmeas-1) - disconnected )