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_complex_math.h>
  17 #include <gsl/gsl_math.h>
  18 #include <gsl/gsl_eigen.h>
  19
  20 using namespace std;
  21
  22 class obs_diagcorr : public o815::obs {
  23
  24 public:
  25   obs_diagcorr(o815 *_O815);
  26
  27 private:
  28   void _start();
  29   void _meas(bool loadedobs, const int& nthmeas);
  30   void _finish();
  31
  32   void corrCompute();
  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);
  35
  36   sim *Sim;
  37
  38   int spatialV;
  39
  40   complex<double> *OM[16+1];
  41
  42   gsl_matrix_complex ***measurements;
  43
  44   gsl_vector_complex **measurements_disco;
  45
  46   static void cdiag (gsl_vector *v, gsl_matrix_complex *m);
  47 };
  48
  49 obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr",
  50                                                     _O815->paraQ->getParaNames() +
  51                                                     "tsep"
  52                                                     "...",
  53                                                     _O815, sizeof(complex<double>) * ( 16*_O815->comargs.lsize[1]/2 + 1*4 ) ) {
  54   for (int ivar = 0; ivar<16+1; ivar++)
  55     OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 ) );
  56
  57   Sim = (sim*)O815->Sim;
  58   spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
  59
  60   measurements = new gsl_matrix_complex**[O815->comargs.nmeas];
  61   for (int imeas = 0; imeas<O815->comargs.nmeas; imeas++) {
  62     measurements[imeas] = new gsl_matrix_complex*[_O815->comargs.lsize[1]/2];
  63     for (int itsep=0; itsep<_O815->comargs.lsize[1]/2; itsep++)
  64       measurements[imeas][itsep] = gsl_matrix_complex_alloc(4,4);
  65   }
  66
  67   measurements_disco = new gsl_vector_complex*[O815->comargs.nmeas];
  68   for (int imeas = 0; imeas<O815->comargs.nmeas; imeas++)
  69     measurements_disco[imeas] = gsl_vector_complex_alloc(4);
  70 }
  71
  72 void obs_diagcorr::_start() {
  73   //  *out << O815->comargs.nmeas << endl;
  74
  75   //*out << "OBS_test: start" << endl;
  76 };
  77
  78 void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
  79   if (!loadedobs)
  80     corrCompute();
  81
  82   for (int icorr=0; icorr<4; icorr++)
  83     for (int jcorr=0; jcorr<4; jcorr++)
  84       for (int itsep=0; itsep<O815->comargs.lsize[1]/2; itsep++) {
  85         gsl_complex tmpc;
  86         tmpc.dat[0] = OM[icorr*4+jcorr][itsep].real();
  87         tmpc.dat[1] = OM[icorr*4+jcorr][itsep].imag();
  88         gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr, tmpc);
  89       }
  90
  91   for (int icorr=0; icorr<4; icorr++) {
  92     gsl_complex tmpc;
  93     tmpc.dat[0] = OM[16][icorr].real();
  94     tmpc.dat[1] = OM[16][icorr].imag();
  95     gsl_vector_complex_set(measurements_disco[nthmeas], icorr, tmpc);
  96   }
  97 };
  98
  99 void obs_diagcorr::cdiag (gsl_vector *v, gsl_matrix_complex *m)
 100 {
 101   gsl_eigen_herm_workspace *wspace = gsl_eigen_herm_alloc(100);
 102   gsl_eigen_herm(m, v, wspace);
 103   gsl_eigen_herm_free(wspace);
 104 }
 105
 106 void obs_diagcorr::_finish() {
 107   gsl_matrix_complex *totalval = gsl_matrix_complex_alloc(4,4);
 108   gsl_matrix_complex *tmpmatrix = gsl_matrix_complex_alloc(4,4);
 109   gsl_vector *tmpvec = gsl_vector_alloc(4);
 110   gsl_vector_complex *tmpvecc = gsl_vector_complex_alloc(4);
 111   gsl_vector *tmpvec2 = gsl_vector_alloc(4);
 112   gsl_vector *jackres = gsl_vector_alloc(4);
 113   gsl_vector *jackerrornorm = gsl_vector_alloc(4);
 114   gsl_vector_complex *totaldisco = gsl_vector_complex_alloc(4);
 115
 116   gsl_vector_complex_set_zero(totaldisco);
 117   for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
 118     gsl_vector_complex_add( totaldisco, measurements_disco[imeas] );
 119
 120   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 121     gsl_matrix_complex_set_zero(totalval);
 122     gsl_vector_set_zero(jackerrornorm);
 123
 124     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
 125       gsl_matrix_complex_add( totalval, measurements[imeas][itsep] );
 126
 127     // compute mean //
 128     gsl_complex tmpc2;
 129     tmpc2.dat[0] = 1.0/O815->comargs.nmeas;
 130     tmpc2.dat[1] = 0;
 131     gsl_matrix_complex_memcpy (tmpmatrix, totalval);
 132     gsl_matrix_complex_scale (tmpmatrix, tmpc2);
 133
 134     for (int icorr=0; icorr<4; icorr++)
 135       for (int jcorr=0; jcorr<4; jcorr++) {
 136         gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(totaldisco, icorr),
 137                                                  gsl_complex_conjugate( gsl_vector_complex_get(totaldisco, jcorr) )
 138                                                  );
 139         discopart = gsl_complex_mul_real(discopart, pow(1.0/O815->comargs.nmeas,2));
 140         gsl_matrix_complex_set( tmpmatrix, icorr, jcorr,
 141                                 gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ),
 142                                                  discopart
 143                                                  )
 144                                 );
 145       }
 146
 147     cdiag(jackres, tmpmatrix);
 148     // END
 149
 150     for (int imeas=0; imeas<O815->comargs.nmeas; imeas++) {
 151       gsl_complex tmpc;
 152       tmpc.dat[0] = 1.0/(O815->comargs.nmeas-1);
 153       tmpc.dat[1] = 0;
 154
 155       gsl_matrix_complex_memcpy (tmpmatrix, totalval);
 156       gsl_matrix_complex_sub (tmpmatrix, measurements[imeas][itsep]);
 157
 158       gsl_matrix_complex_scale ( tmpmatrix, tmpc );
 159
 160       gsl_vector_complex_memcpy( tmpvecc, totaldisco );
 161       gsl_vector_complex_sub( tmpvecc, measurements_disco[imeas] );
 162
 163       for (int icorr=0; icorr<4; icorr++)
 164         for (int jcorr=0; jcorr<4; jcorr++) {
 165           gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(tmpvecc, icorr),
 166                                                    gsl_complex_conjugate( gsl_vector_complex_get(tmpvecc, jcorr) )
 167                                                    );
 168           discopart = gsl_complex_mul_real(discopart, pow(1.0/(O815->comargs.nmeas-1),2));
 169           gsl_matrix_complex_set( tmpmatrix, icorr, jcorr,
 170                                   gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ),
 171                                                    discopart
 172                                                    )
 173                                   );
 174         }
 175
 176       cdiag(tmpvec, tmpmatrix);
 177
 178       gsl_vector_sub(tmpvec, jackres);
 179       gsl_vector_memcpy(tmpvec2, tmpvec);
 180       gsl_vector_mul(tmpvec, tmpvec2);
 181
 182       gsl_vector_add(jackerrornorm, tmpvec);
 183     }
 184
 185     gsl_vector_scale( jackerrornorm,
 186                       (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas );
 187
 188     *out << O815->paraQ->getParaVals();
 189     *out << "\t" << itsep;
 190     for (int iev = 0; iev < 4; iev++) {
 191       *out << "\t" << gsl_vector_get(jackres,iev)
 192            << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev));
 193     }
 194     *out << endl;
 195   }
 196
 197   gsl_matrix_complex_free(totalval);
 198 };
 199
 200 void obs_diagcorr::corrCompute()
 201 {
 202   complex<double> phislice[4][O815->comargs.lsize[1]];
 203
 204   for (int icorr=0; icorr<4; icorr++)
 205     OM[16][icorr] = 0;
 206
 207   for (int it = 0; it < O815->comargs.lsize[1]; it++) {
 208     for (int icorr=0; icorr<4; icorr++)
 209       phislice[icorr][it] = 0;
 210
 211     for (int ix = 0; ix < spatialV; ix++) {
 212       phislice[0][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 213       phislice[1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ] * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 214       phislice[2][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]);
 215       phislice[3][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 216     }
 217
 218     for (int icorr=0; icorr<4; icorr++) {
 219       phislice[icorr][it] /= spatialV;
 220       OM[16][icorr] += phislice[icorr][it];
 221     }
 222   }
 223
 224   for (int icorr = 0; icorr < 4; icorr++) {
 225     for (int jcorr = 0; jcorr < 4; jcorr++) {
 226       for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 227         OM[icorr*4+jcorr][itsep] = 0;
 228
 229         for (int it = 0; it < O815->comargs.lsize[1]; it++)
 230           OM[icorr*4+jcorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[jcorr][it] );
 231
 232         OM[icorr*4+jcorr][itsep] /= O815->comargs.lsize[1];
 233       }
 234     }
 235
 236     OM[16][icorr] /= O815->comargs.lsize[1];
 237   }
 238 }
 239
 240 #endif