...

diff --git a/u1casc-ordinary/obs_diagcorr.hpp b/u1casc-ordinary/obs_diagcorr.hpp
new file mode 100644 (file)
index 0000000..e5701fd
--- /dev/null
+++ b/u1casc-ordinary/obs_diagcorr.hpp
@@ -0,0 +1,159 @@
+#ifndef OBS_DIAGCORR_HPP
+#define OBS_DIAGCORR_HPP
+
+#include "latlib/o815/o815.h"
+
+#include "latlib/writeout.h"
+
+#include "latlib/obstat.hpp"
+
+#include <iostream>
+#include <complex>
+
+#include <cmath>
+
+#include <gsl/gsl_complex.h>
+#include <gsl/gsl_math.h>
+#include <gsl/gsl_eigen.h>
+
+using namespace std;
+
+class obs_diagcorr : public o815::obs {
+
+public:
+  obs_diagcorr(o815 *_O815);
+  
+private:
+  void _start();
+  void _meas(bool loadedobs, const int& nthmeas);
+  void _finish();
+
+  void corrCompute();
+  static complex<double> effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para);
+  static void preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para);
+
+  sim *Sim;
+  obstat< complex<double>, complex<double>  > oC[16];
+
+  int spatialV;
+
+  complex<double> *OM[16];
+};
+
+obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr", 
+                                           _O815->paraQ->getParaNames() + 
+                                           "tsep"
+                                           "...", 
+                                           _O815, 16*sizeof(complex<double>)*(_O815->comargs.lsize[1]/2) ) {
+  for (int ivar = 0; ivar<16; ivar++)
+    OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 ) );
+  
+  Sim = (sim*)O815->Sim;
+  spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
+}
+
+void obs_diagcorr::_start() {
+  //*out << "OBS_test: start" << endl;
+};
+
+void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
+  if (!loadedobs)
+    corrCompute();
+
+  for (int icorr=0; icorr<16; icorr++)
+    oC[icorr].addMeas( OM[icorr], O815->comargs.lsize[1]/2 );
+};
+
+void obs_diagcorr::_finish() {
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+    gsl_complex matdata[16];
+
+    gsl_matrix_complex *m = gsl_matrix_complex_alloc(4,4);
+    
+    for (int icorr = 0; icorr < 4; icorr++) 
+      for (int jcorr = 0; jcorr < 4; jcorr++) {
+       int compid_corr = oC[icorr*4+jcorr].computeEasy(itsep);
+       gsl_complex tmpc;
+       
+       tmpc.dat[0] = oC[icorr*4+jcorr].getMean(compid_corr).real();
+       tmpc.dat[1] = oC[icorr*4+jcorr].getMean(compid_corr).imag();
+       
+       gsl_matrix_complex_set( m, icorr, jcorr, tmpc);
+      }
+
+    gsl_vector *evvec = gsl_vector_alloc(4);
+
+    gsl_eigen_herm_workspace *wspace = gsl_eigen_herm_alloc(100);
+
+    gsl_eigen_herm(m, evvec, wspace);
+
+    gsl_eigen_herm_free(wspace);
+    gsl_matrix_complex_free(m);
+    
+    *out << O815->paraQ->getParaVals();
+    *out << "\t" << itsep;
+    for (int iev = 0; iev < 4; iev++) {
+      *out << "\t" << gsl_vector_get(evvec,iev);
+    }
+    *out << endl;
+
+    gsl_vector_free(evvec);
+  }
+  
+  
+  // for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+  //   complex
+  // }
+
+    /*
+      for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+    *out << O815->paraQ->getParaVals();
+    *out << "\t" << itsep;
+    for (int icorr = 0; icorr < 16; icorr++) {
+      *out << "\t" << real( oC[icorr].getMean(compid_corr[itsep][icorr]) );
+    }
+    *out << endl;
+    }
+    */
+
+  for (int icorr = 0; icorr < 16; icorr++) {
+    oC[icorr].reset();
+  }
+};
+
+void obs_diagcorr::corrCompute()
+{
+  complex<double> phislice[4][O815->comargs.lsize[1]];
+
+  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 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 ];
+      phislice[1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ] * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
+      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++)
+    for (int jcorr = 0; jcorr < 4; jcorr++) {
+      for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+       OM[icorr*4+jcorr][itsep] = 0;
+       
+       for (int it = 0; it < O815->comargs.lsize[1]; it++)
+         OM[icorr*4+jcorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[jcorr][it] );
+      
+       OM[icorr*4+jcorr][itsep] /= O815->comargs.lsize[1];
+      }
+    }
+}
+
+#endif