X-Git-Url: http://git.treefish.org/~alex/phys/u1casc.git/blobdiff_plain/33f700a8121cd3f356902483d7af9a9a253037a0..607713e7bb60eb68ad4285a805f10943a7cec2ff:/u1casc-ordinary/obs_diagcorr.hpp diff --git a/u1casc-ordinary/obs_diagcorr.hpp b/u1casc-ordinary/obs_diagcorr.hpp index 43d2132..a8d88b5 100644 --- a/u1casc-ordinary/obs_diagcorr.hpp +++ b/u1casc-ordinary/obs_diagcorr.hpp @@ -16,6 +16,9 @@ #include #include #include +#include +#include +#include using namespace std; @@ -43,7 +46,7 @@ private: gsl_vector_complex **measurements_disco; - static void cdiag (gsl_vector *v, gsl_matrix_complex *m); + static void cdiag (gsl_vector *v, gsl_matrix_complex *evec, gsl_matrix_complex *m); }; obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr", @@ -82,37 +85,40 @@ 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_complex tmpc; - tmpc.dat[0] = OM[16][icorr].real(); - tmpc.dat[1] = OM[16][icorr].imag(); - gsl_vector_complex_set(measurements_disco[nthmeas], icorr, tmpc); + 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) +void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *evec, gsl_matrix_complex *m) { - gsl_eigen_herm_workspace *wspace = gsl_eigen_herm_alloc(100); - gsl_eigen_herm(m, v, wspace); - gsl_eigen_herm_free(wspace); + gsl_eigen_hermv_workspace *wspace = gsl_eigen_hermv_alloc(100); + + gsl_eigen_hermv(m, v, evec, wspace); + gsl_eigen_hermv_sort(v, evec, GSL_EIGEN_SORT_VAL_ASC); + + gsl_eigen_hermv_free(wspace); + gsl_matrix_complex_free(evec); } 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_matrix_complex *tmpmatrix2 = 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_vector_complex *totaldisco = gsl_vector_complex_alloc(4); - + gsl_matrix_complex *evecres = gsl_matrix_complex_alloc(4,4); + gsl_matrix_complex *evecerrornorm = gsl_matrix_complex_alloc(4,4); + gsl_vector_complex_set_zero(totaldisco); for (int imeas=0; imeascomargs.nmeas; imeas++) gsl_vector_complex_add( totaldisco, measurements_disco[imeas] ); @@ -120,17 +126,13 @@ void obs_diagcorr::_finish() { for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) { gsl_matrix_complex_set_zero(totalval); gsl_vector_set_zero(jackerrornorm); + gsl_matrix_complex_set_zero(evecerrornorm); 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), @@ -143,23 +145,14 @@ void obs_diagcorr::_finish() { ) ); } - - cdiag(jackres, tmpmatrix); - // END + cdiag(jackres, evecres, tmpmatrix); 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), @@ -171,30 +164,49 @@ void obs_diagcorr::_finish() { discopart ) ); - } + } + cdiag(tmpvec, tmpmatrix2, tmpmatrix); - cdiag(tmpvec, tmpmatrix); - gsl_vector_sub(tmpvec, jackres); gsl_vector_memcpy(tmpvec2, tmpvec); gsl_vector_mul(tmpvec, tmpvec2); - gsl_vector_add(jackerrornorm, tmpvec); - } - + + gsl_matrix_complex_sub(tmpmatrix2, evecres); + for (int i=0; i<4; i++) + for (int j=0; j<4; j++) { + gsl_complex tmperr = + gsl_complex_rect( pow(GSL_REAL( gsl_matrix_complex_get( tmpmatrix2, i, j ) ),2), + pow(GSL_IMAG( gsl_matrix_complex_get( tmpmatrix2, i, j ) ),2)); + gsl_matrix_complex_set(tmpmatrix, i, j, tmperr); + } + gsl_matrix_complex_add( evecerrornorm, tmpmatrix ); + } gsl_vector_scale( jackerrornorm, (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas ); - + gsl_matrix_complex_scale( evecerrornorm, gsl_complex_rect( (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas, 0.0 ) ); + *out << O815->paraQ->getParaVals(); *out << "\t" << itsep; for (int iev = 0; iev < 4; iev++) { *out << "\t" << gsl_vector_get(jackres,iev) << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev)); + for (int iinter=0; iinter<4; iinter++) + *out << "\t" << GSL_REAL(gsl_matrix_complex_get(evecres, iinter, iev)) + << "\t" << sqrt( GSL_REAL(gsl_matrix_complex_get(evecerrornorm, iinter, iev)) ); } *out << endl; } gsl_matrix_complex_free(totalval); + gsl_matrix_complex_free(tmpmatrix); + gsl_matrix_complex_free(evecres); + 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()