X-Git-Url: http://git.treefish.org/~alex/phys/u1casc.git/blobdiff_plain/ea278d0dfaf581ddaf7b1b16af83dc1b9b62b754..65fa392bb1e02b52b05bea7b7c10efc6f421830a:/u1casc-ordinary/obs_diagcorr.hpp diff --git a/u1casc-ordinary/obs_diagcorr.hpp b/u1casc-ordinary/obs_diagcorr.hpp index 3ba3960..57b835b 100644 --- a/u1casc-ordinary/obs_diagcorr.hpp +++ b/u1casc-ordinary/obs_diagcorr.hpp @@ -13,6 +13,7 @@ #include #include +#include #include #include @@ -33,23 +34,24 @@ private: static void preEffMass(vector< vector < complex > > *allVals, vector < complex > *preCalculated, void *para); sim *Sim; - obstat< complex, complex > oC[16]; int spatialV; - complex *OM[16]; + complex *OM[16+1]; gsl_matrix_complex ***measurements; + gsl_vector_complex **measurements_disco; + static void cdiag (gsl_vector *v, gsl_matrix_complex *m); }; obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr", - _O815->paraQ->getParaNames() + - "tsep" - "...", - _O815, 16*sizeof(complex)*(_O815->comargs.lsize[1]/2) ) { - for (int ivar = 0; ivar<16; ivar++) + _O815->paraQ->getParaNames() + + "tsep" + "...", + _O815, sizeof(complex) * ( 16*_O815->comargs.lsize[1]/2 + 1*4 ) ) { + for (int ivar = 0; ivar<16+1; ivar++) OM[ivar] = (complex*)( obsMem + ivar*sizeof(complex)*( _O815->comargs.lsize[1]/2 ) ); Sim = (sim*)O815->Sim; @@ -61,6 +63,10 @@ obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr", for (int itsep=0; itsep<_O815->comargs.lsize[1]/2; itsep++) measurements[imeas][itsep] = gsl_matrix_complex_alloc(4,4); } + + measurements_disco = new gsl_vector_complex*[O815->comargs.nmeas]; + for (int imeas = 0; imeascomargs.nmeas; imeas++) + measurements_disco[imeas] = gsl_vector_complex_alloc(4); } void obs_diagcorr::_start() { @@ -76,11 +82,14 @@ void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) { for (int icorr=0; icorr<4; icorr++) for (int jcorr=0; jcorr<4; jcorr++) for (int itsep=0; itsepcomargs.lsize[1]/2; itsep++) { - gsl_complex tmpc; - tmpc.dat[0] = OM[icorr*4+jcorr][itsep].real(); - tmpc.dat[1] = OM[icorr*4+jcorr][itsep].imag(); - gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr, tmpc); + gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr, + gsl_complex_rect( OM[icorr*4+jcorr][itsep].real(), OM[icorr*4+jcorr][itsep].imag() )); } + + for (int icorr=0; icorr<4; icorr++) { + gsl_vector_complex_set(measurements_disco[nthmeas], icorr, + gsl_complex_rect( OM[16][icorr].real(), OM[16][icorr].imag() ) ); + } }; void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *m) @@ -94,50 +103,63 @@ void obs_diagcorr::_finish() { gsl_matrix_complex *totalval = gsl_matrix_complex_alloc(4,4); gsl_matrix_complex *tmpmatrix = gsl_matrix_complex_alloc(4,4); gsl_vector *tmpvec = gsl_vector_alloc(4); + gsl_vector_complex *tmpvecc = gsl_vector_complex_alloc(4); gsl_vector *tmpvec2 = gsl_vector_alloc(4); gsl_vector *jackres = gsl_vector_alloc(4); gsl_vector *jackerrornorm = gsl_vector_alloc(4); - //gsl_matrix_complex *manymeans[O815->comargs.nmeas]; - - /* - for (int imeas=0; imeascomargs.lsize[1]/2; imeas++) - manymeans[imeas] = gsl_matrix_complex_alloc(4,4); - */ + gsl_vector_complex *totaldisco = gsl_vector_complex_alloc(4); + gsl_vector_complex_set_zero(totaldisco); + for (int imeas=0; imeascomargs.nmeas; imeas++) + gsl_vector_complex_add( totaldisco, measurements_disco[imeas] ); + for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) { gsl_matrix_complex_set_zero(totalval); gsl_vector_set_zero(jackerrornorm); for (int imeas=0; imeascomargs.nmeas; imeas++) gsl_matrix_complex_add( totalval, measurements[imeas][itsep] ); - - // compute mean // - gsl_complex tmpc2; - tmpc2.dat[0] = 1.0/O815->comargs.nmeas; - tmpc2.dat[1] = 0; + gsl_matrix_complex_memcpy (tmpmatrix, totalval); - gsl_matrix_complex_scale (tmpmatrix, tmpc2); + gsl_matrix_complex_scale (tmpmatrix, gsl_complex_rect(1.0/O815->comargs.nmeas, 0.0) ); + for (int icorr=0; icorr<4; icorr++) + for (int jcorr=0; jcorr<4; jcorr++) { + gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(totaldisco, icorr), + gsl_complex_conjugate( gsl_vector_complex_get(totaldisco, jcorr) ) + ); + discopart = gsl_complex_mul_real(discopart, pow(1.0/O815->comargs.nmeas,2)); + gsl_matrix_complex_set( tmpmatrix, icorr, jcorr, + gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ), + discopart + ) + ); + } cdiag(jackres, tmpmatrix); - // END for (int imeas=0; imeascomargs.nmeas; imeas++) { - gsl_complex tmpc; - tmpc.dat[0] = 1.0/(O815->comargs.nmeas-1); - tmpc.dat[1] = 0; - gsl_matrix_complex_memcpy (tmpmatrix, totalval); gsl_matrix_complex_sub (tmpmatrix, measurements[imeas][itsep]); - gsl_matrix_complex_scale ( tmpmatrix, tmpc ); - + gsl_matrix_complex_scale ( tmpmatrix, gsl_complex_rect(1.0/(O815->comargs.nmeas-1), 0.0) ); + gsl_vector_complex_memcpy( tmpvecc, totaldisco ); + gsl_vector_complex_sub( tmpvecc, measurements_disco[imeas] ); + for (int icorr=0; icorr<4; icorr++) + for (int jcorr=0; jcorr<4; jcorr++) { + gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(tmpvecc, icorr), + gsl_complex_conjugate( gsl_vector_complex_get(tmpvecc, jcorr) ) + ); + discopart = gsl_complex_mul_real(discopart, pow(1.0/(O815->comargs.nmeas-1),2)); + gsl_matrix_complex_set( tmpmatrix, icorr, jcorr, + gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ), + discopart + ) + ); + } cdiag(tmpvec, tmpmatrix); - gsl_vector_sub(tmpvec, jackres); gsl_vector_memcpy(tmpvec2, tmpvec); - gsl_vector_mul(tmpvec, tmpvec2); //!!!!!! - + gsl_vector_mul(tmpvec, tmpvec2); gsl_vector_add(jackerrornorm, tmpvec); - } - + } gsl_vector_scale( jackerrornorm, (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas ); @@ -148,24 +170,28 @@ void obs_diagcorr::_finish() { << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev)); } *out << endl; - - //for (int imeas=0; imeascomargs.nmeas; imeas++) - //gsl_matrix_ - // gsl_matrix_complex_memcpy( manymeans[imeas], precalc ); } gsl_matrix_complex_free(totalval); + gsl_matrix_complex_free(tmpmatrix); + gsl_vector_free(tmpvec); + gsl_vector_complex_free(tmpvecc); + gsl_vector_free(tmpvec2); + gsl_vector_free(jackres); + gsl_vector_free(jackerrornorm); + gsl_vector_complex_free(totaldisco); }; void obs_diagcorr::corrCompute() { complex phislice[4][O815->comargs.lsize[1]]; + for (int icorr=0; icorr<4; icorr++) + OM[16][icorr] = 0; + for (int it = 0; it < O815->comargs.lsize[1]; it++) { - phislice[0][it] = 0; - phislice[1][it] = 0; - phislice[2][it] = 0; - phislice[3][it] = 0; + for (int icorr=0; icorr<4; icorr++) + phislice[icorr][it] = 0; for (int ix = 0; ix < spatialV; ix++) { phislice[0][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]; @@ -173,13 +199,14 @@ void obs_diagcorr::corrCompute() phislice[2][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]); phislice[3][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]; } - phislice[0][it] /= spatialV; - phislice[1][it] /= spatialV; - phislice[2][it] /= spatialV; - phislice[3][it] /= spatialV; + + for (int icorr=0; icorr<4; icorr++) { + phislice[icorr][it] /= spatialV; + OM[16][icorr] += phislice[icorr][it]; + } } - for (int icorr = 0; icorr < 4; icorr++) + for (int icorr = 0; icorr < 4; icorr++) { for (int jcorr = 0; jcorr < 4; jcorr++) { for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) { OM[icorr*4+jcorr][itsep] = 0; @@ -190,6 +217,9 @@ void obs_diagcorr::corrCompute() OM[icorr*4+jcorr][itsep] /= O815->comargs.lsize[1]; } } + + OM[16][icorr] /= O815->comargs.lsize[1]; + } } #endif