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
  43 obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr",
  44                                             _O815->paraQ->getParaNames() +
  45                                             "tsep"
  46                                             "...",
  47                                             _O815, 16*sizeof(complex<double>)*(_O815->comargs.lsize[1]/2) ) {
  48   for (int ivar = 0; ivar<16; ivar++)
  49     OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 ) );
  50
  51   Sim = (sim*)O815->Sim;
  52   spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
  53 }
  54
  55 void obs_diagcorr::_start() {
  56   //*out << "OBS_test: start" << endl;
  57 };
  58
  59 void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
  60   if (!loadedobs)
  61     corrCompute();
  62
  63   for (int icorr=0; icorr<16; icorr++)
  64     oC[icorr].addMeas( OM[icorr], O815->comargs.lsize[1]/2 );
  65 };
  66
  67 void obs_diagcorr::_finish() {
  68   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
  69     gsl_complex matdata[16];
  70
  71     gsl_matrix_complex *m = gsl_matrix_complex_alloc(4,4);
  72
  73     for (int icorr = 0; icorr < 4; icorr++)
  74       for (int jcorr = 0; jcorr < 4; jcorr++) {
  75         int compid_corr = oC[icorr*4+jcorr].computeEasy(itsep);
  76         gsl_complex tmpc;
  77
  78         tmpc.dat[0] = oC[icorr*4+jcorr].getMean(compid_corr).real();
  79         tmpc.dat[1] = oC[icorr*4+jcorr].getMean(compid_corr).imag();
  80
  81         gsl_matrix_complex_set( m, icorr, jcorr, tmpc);
  82       }
  83
  84     gsl_vector *evvec = gsl_vector_alloc(4);
  85
  86     gsl_eigen_herm_workspace *wspace = gsl_eigen_herm_alloc(100);
  87
  88     gsl_eigen_herm(m, evvec, wspace);
  89
  90     gsl_eigen_herm_free(wspace);
  91     gsl_matrix_complex_free(m);
  92
  93     *out << O815->paraQ->getParaVals();
  94     *out << "\t" << itsep;
  95     for (int iev = 0; iev < 4; iev++) {
  96       *out << "\t" << gsl_vector_get(evvec,iev);
  97     }
  98     *out << endl;
  99
 100     gsl_vector_free(evvec);
 101   }
 102
 103
 104   // for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 105   //   complex
 106   // }
 107
 108     /*
 109       for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 110     *out << O815->paraQ->getParaVals();
 111     *out << "\t" << itsep;
 112     for (int icorr = 0; icorr < 16; icorr++) {
 113       *out << "\t" << real( oC[icorr].getMean(compid_corr[itsep][icorr]) );
 114     }
 115     *out << endl;
 116     }
 117     */
 118
 119   for (int icorr = 0; icorr < 16; icorr++) {
 120     oC[icorr].reset();
 121   }
 122 };
 123
 124 void obs_diagcorr::corrCompute()
 125 {
 126   complex<double> phislice[4][O815->comargs.lsize[1]];
 127
 128   for (int it = 0; it < O815->comargs.lsize[1]; it++) {
 129     phislice[0][it] = 0;
 130     phislice[1][it] = 0;
 131     phislice[2][it] = 0;
 132     phislice[3][it] = 0;
 133
 134     for (int ix = 0; ix < spatialV; ix++) {
 135       phislice[0][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 136       phislice[1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ] * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 137       phislice[2][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]);
 138       phislice[3][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 139     }
 140     phislice[0][it] /= spatialV;
 141     phislice[1][it] /= spatialV;
 142     phislice[2][it] /= spatialV;
 143     phislice[3][it] /= spatialV;
 144   }
 145
 146   for (int icorr = 0; icorr < 4; icorr++)
 147     for (int jcorr = 0; jcorr < 4; jcorr++) {
 148       for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 149         OM[icorr*4+jcorr][itsep] = 0;
 150
 151         for (int it = 0; it < O815->comargs.lsize[1]; it++)
 152           OM[icorr*4+jcorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[jcorr][it] );
 153
 154         OM[icorr*4+jcorr][itsep] /= O815->comargs.lsize[1];
 155       }
 156     }
 157 }
 158
 159 #endif