...

[phys/u1casc.git] / u1casc-ordinary / obs_diagcorr.hpp

diff --git a/u1casc-ordinary/obs_diagcorr.hpp b/u1casc-ordinary/obs_diagcorr.hpp
index e5701fdff714e0829205780482208d7b2b01f3f2..3ba396097f1e457b50705e91cb4557c110791836 100644 (file)
--- a/u1casc-ordinary/obs_diagcorr.hpp
+++ b/u1casc-ordinary/obs_diagcorr.hpp
@@ -38,6 +38,10 @@ private:
   int spatialV;
 
   complex<double> *OM[16];
+
+  gsl_matrix_complex ***measurements;
+
+  static void cdiag (gsl_vector *v, gsl_matrix_complex *m);
 };
 
 obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr", 
@@ -50,9 +54,18 @@ obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr",
   
   Sim = (sim*)O815->Sim;
   spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
+
+  measurements = new gsl_matrix_complex**[O815->comargs.nmeas];
+  for (int imeas = 0; imeas<O815->comargs.nmeas; imeas++) {
+    measurements[imeas] = new gsl_matrix_complex*[_O815->comargs.lsize[1]/2];
+    for (int itsep=0; itsep<_O815->comargs.lsize[1]/2; itsep++)
+      measurements[imeas][itsep] = gsl_matrix_complex_alloc(4,4);
+  }
 }
 
 void obs_diagcorr::_start() {
+  //  *out << O815->comargs.nmeas << endl;
+  
   //*out << "OBS_test: start" << endl;
 };
 
@@ -60,65 +73,88 @@ 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 );
+  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++) {
+       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);
+      }
 };
 
-void obs_diagcorr::_finish() {
-  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
-    gsl_complex matdata[16];
+void obs_diagcorr::cdiag (gsl_vector *v, 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_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);
-      }
+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_vector *tmpvec = gsl_vector_alloc(4);
+  gsl_vector *tmpvec2 = gsl_vector_alloc(4);
+  gsl_vector *jackres = gsl_vector_alloc(4);
+  gsl_vector *jackerrornorm = gsl_vector_alloc(4);
+  //gsl_matrix_complex *manymeans[O815->comargs.nmeas];
+
+  /*
+  for (int imeas=0; imeas<O815->comargs.lsize[1]/2; imeas++)
+    manymeans[imeas] = gsl_matrix_complex_alloc(4,4);
+  */
 
-    gsl_vector *evvec = gsl_vector_alloc(4);
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+    gsl_matrix_complex_set_zero(totalval);
+    gsl_vector_set_zero(jackerrornorm);
+
+    for (int imeas=0; imeas<O815->comargs.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);
+    cdiag(jackres, tmpmatrix);
+    // END
+
+    for (int imeas=0; imeas<O815->comargs.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_eigen_herm_workspace *wspace = gsl_eigen_herm_alloc(100);
+      cdiag(tmpvec, tmpmatrix);
 
-    gsl_eigen_herm(m, evvec, wspace);
+      gsl_vector_sub(tmpvec, jackres);
+      gsl_vector_memcpy(tmpvec2, tmpvec);
+      gsl_vector_mul(tmpvec, tmpvec2); //!!!!!!
 
-    gsl_eigen_herm_free(wspace);
-    gsl_matrix_complex_free(m);
+      gsl_vector_add(jackerrornorm, tmpvec);
+    }
     
+    gsl_vector_scale( jackerrornorm,
+                     (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas );
+
     *out << O815->paraQ->getParaVals();
     *out << "\t" << itsep;
     for (int iev = 0; iev < 4; iev++) {
-      *out << "\t" << gsl_vector_get(evvec,iev);
+      *out << "\t" << gsl_vector_get(jackres,iev)
+          << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev));
     }
     *out << endl;
-
-    gsl_vector_free(evvec);
+    
+    //for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
+      //gsl_matrix_
+      //  gsl_matrix_complex_memcpy( manymeans[imeas], precalc );
   }
-  
-  
-  // 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();
-  }
+  gsl_matrix_complex_free(totalval);
 };
 
 void obs_diagcorr::corrCompute()