u1casc-ordinary/obs_diagcorr.hpp

   1 #ifndef OBS_DIAGCORR_HPP
   2 #define OBS_DIAGCORR_HPP
   3
   4 #include "latlib/o815/o815.h"
   5
   6 #include "latlib/writeout.h"
   7
   8 #include "latlib/obstat.hpp"
   9
  10 #include <iostream>
  11 #include <complex>
  12
  13 #include <cmath>
  14
  15 #include <gsl/gsl_complex.h>
  16 #include <gsl/gsl_math.h>
  17 #include <gsl/gsl_eigen.h>
  18
  19 using namespace std;
  20
  21 class obs_diagcorr : public o815::obs {
  22
  23 public:
  24   obs_diagcorr(o815 *_O815);
  25
  26 private:
  27   void _start();
  28   void _meas(bool loadedobs, const int& nthmeas);
  29   void _finish();
  30
  31   void corrCompute();
  32   static complex<double> effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para);
  33   static void preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para);
  34
  35   sim *Sim;
  36   obstat< complex<double>, complex<double>  > oC[16];
  37
  38   int spatialV;
  39
  40   complex<double> *OM[16];
  41
  42   gsl_matrix_complex ***measurements;
  43
  44   static void cdiag (gsl_vector *v, gsl_matrix_complex *m);
  45 };
  46
  47 obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr",
  48                                             _O815->paraQ->getParaNames() +
  49                                             "tsep"
  50                                             "...",
  51                                             _O815, 16*sizeof(complex<double>)*(_O815->comargs.lsize[1]/2) ) {
  52   for (int ivar = 0; ivar<16; ivar++)
  53     OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 ) );
  54
  55   Sim = (sim*)O815->Sim;
  56   spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
  57
  58   measurements = new gsl_matrix_complex**[O815->comargs.nmeas];
  59   for (int imeas = 0; imeas<O815->comargs.nmeas; imeas++) {
  60     measurements[imeas] = new gsl_matrix_complex*[_O815->comargs.lsize[1]/2];
  61     for (int itsep=0; itsep<_O815->comargs.lsize[1]/2; itsep++)
  62       measurements[imeas][itsep] = gsl_matrix_complex_alloc(4,4);
  63   }
  64 }
  65
  66 void obs_diagcorr::_start() {
  67   //  *out << O815->comargs.nmeas << endl;
  68
  69   //*out << "OBS_test: start" << endl;
  70 };
  71
  72 void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
  73   if (!loadedobs)
  74     corrCompute();
  75
  76   for (int icorr=0; icorr<4; icorr++)
  77     for (int jcorr=0; jcorr<4; jcorr++)
  78       for (int itsep=0; itsep<O815->comargs.lsize[1]/2; itsep++) {
  79         gsl_complex tmpc;
  80         tmpc.dat[0] = OM[icorr*4+jcorr][itsep].real();
  81         tmpc.dat[1] = OM[icorr*4+jcorr][itsep].imag();
  82         gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr, tmpc);
  83       }
  84 };
  85
  86 void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *m)
  87 {
  88   gsl_eigen_herm_workspace *wspace = gsl_eigen_herm_alloc(100);
  89   gsl_eigen_herm(m, v, wspace);
  90   gsl_eigen_herm_free(wspace);
  91 }
  92
  93 void obs_diagcorr::_finish() {
  94   gsl_matrix_complex *totalval = gsl_matrix_complex_alloc(4,4);
  95   gsl_matrix_complex *tmpmatrix = gsl_matrix_complex_alloc(4,4);
  96   gsl_vector *tmpvec = gsl_vector_alloc(4);
  97   gsl_vector *tmpvec2 = gsl_vector_alloc(4);
  98   gsl_vector *jackres = gsl_vector_alloc(4);
  99   gsl_vector *jackerrornorm = gsl_vector_alloc(4);
 100   //gsl_matrix_complex *manymeans[O815->comargs.nmeas];
 101
 102   /*
 103   for (int imeas=0; imeas<O815->comargs.lsize[1]/2; imeas++)
 104     manymeans[imeas] = gsl_matrix_complex_alloc(4,4);
 105   */
 106
 107   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 108     gsl_matrix_complex_set_zero(totalval);
 109     gsl_vector_set_zero(jackerrornorm);
 110
 111     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
 112       gsl_matrix_complex_add( totalval, measurements[imeas][itsep] );
 113
 114     // compute mean //
 115     gsl_complex tmpc2;
 116     tmpc2.dat[0] = 1.0/O815->comargs.nmeas;
 117     tmpc2.dat[1] = 0;
 118     gsl_matrix_complex_memcpy (tmpmatrix, totalval);
 119     gsl_matrix_complex_scale (tmpmatrix, tmpc2);
 120     cdiag(jackres, tmpmatrix);
 121     // END
 122
 123     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++) {
 124       gsl_complex tmpc;
 125       tmpc.dat[0] = 1.0/(O815->comargs.nmeas-1);
 126       tmpc.dat[1] = 0;
 127
 128       gsl_matrix_complex_memcpy (tmpmatrix, totalval);
 129       gsl_matrix_complex_sub (tmpmatrix, measurements[imeas][itsep]);
 130       gsl_matrix_complex_scale ( tmpmatrix, tmpc );
 131
 132       cdiag(tmpvec, tmpmatrix);
 133
 134       gsl_vector_sub(tmpvec, jackres);
 135       gsl_vector_memcpy(tmpvec2, tmpvec);
 136       gsl_vector_mul(tmpvec, tmpvec2); //!!!!!!
 137
 138       gsl_vector_add(jackerrornorm, tmpvec);
 139     }
 140
 141     gsl_vector_scale( jackerrornorm,
 142                       (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas );
 143
 144     *out << O815->paraQ->getParaVals();
 145     *out << "\t" << itsep;
 146     for (int iev = 0; iev < 4; iev++) {
 147       *out << "\t" << gsl_vector_get(jackres,iev)
 148            << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev));
 149     }
 150     *out << endl;
 151
 152     //for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
 153       //gsl_matrix_
 154       //  gsl_matrix_complex_memcpy( manymeans[imeas], precalc );
 155   }
 156
 157   gsl_matrix_complex_free(totalval);
 158 };
 159
 160 void obs_diagcorr::corrCompute()
 161 {
 162   complex<double> phislice[4][O815->comargs.lsize[1]];
 163
 164   for (int it = 0; it < O815->comargs.lsize[1]; it++) {
 165     phislice[0][it] = 0;
 166     phislice[1][it] = 0;
 167     phislice[2][it] = 0;
 168     phislice[3][it] = 0;
 169
 170     for (int ix = 0; ix < spatialV; ix++) {
 171       phislice[0][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 172       phislice[1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ] * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 173       phislice[2][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]);
 174       phislice[3][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 175     }
 176     phislice[0][it] /= spatialV;
 177     phislice[1][it] /= spatialV;
 178     phislice[2][it] /= spatialV;
 179     phislice[3][it] /= spatialV;
 180   }
 181
 182   for (int icorr = 0; icorr < 4; icorr++)
 183     for (int jcorr = 0; jcorr < 4; jcorr++) {
 184       for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 185         OM[icorr*4+jcorr][itsep] = 0;
 186
 187         for (int it = 0; it < O815->comargs.lsize[1]; it++)
 188           OM[icorr*4+jcorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[jcorr][it] );
 189
 190         OM[icorr*4+jcorr][itsep] /= O815->comargs.lsize[1];
 191       }
 192     }
 193 }
 194
 195 #endif