X-Git-Url: http://git.treefish.org/~alex/phys/u1casc.git/blobdiff_plain/ea278d0dfaf581ddaf7b1b16af83dc1b9b62b754..a47125a5c94fc75125547d64c693e7384ed9b70d:/u1casc-ordinary/obs_diagcorr.hpp

diff --git a/u1casc-ordinary/obs_diagcorr.hpp b/u1casc-ordinary/obs_diagcorr.hpp
index 3ba3960..f51df7c 100644
--- a/u1casc-ordinary/obs_diagcorr.hpp
+++ b/u1casc-ordinary/obs_diagcorr.hpp
@@ -13,8 +13,12 @@
 #include <cmath>
 
 #include <gsl/gsl_complex.h>
+#include <gsl/gsl_complex_math.h>
 #include <gsl/gsl_math.h>
 #include <gsl/gsl_eigen.h>
+#include <gsl/gsl_linalg.h>
+#include <gsl/gsl_cblas.h>
+#include <gsl/gsl_blas.h>
 
 using namespace std;
 
@@ -33,23 +37,26 @@ private:
   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];
+  complex<double> *OM[16+1];
 
   gsl_matrix_complex ***measurements;
 
+  gsl_vector_complex **measurements_disco;
+  
   static void cdiag (gsl_vector *v, gsl_matrix_complex *m);
+
+  static void invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m);
 };
 
 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++)
+						    _O815->paraQ->getParaNames() + 
+						    "tsep"
+						    "...", 
+						    _O815, sizeof(complex<double>) * ( 16*_O815->comargs.lsize[1]/2 + 1*4 ) ) {
+  for (int ivar = 0; ivar<16+1; ivar++)
     OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 ) );
   
   Sim = (sim*)O815->Sim;
@@ -61,6 +68,10 @@ obs_diagcorr::obs_diagcorr(o815 *_O815) : o815::obs("diagcorr",
     for (int itsep=0; itsep<_O815->comargs.lsize[1]/2; itsep++)
       measurements[imeas][itsep] = gsl_matrix_complex_alloc(4,4);
   }
+
+  measurements_disco = new gsl_vector_complex*[O815->comargs.nmeas];
+  for (int imeas = 0; imeas<O815->comargs.nmeas; imeas++)
+    measurements_disco[imeas] = gsl_vector_complex_alloc(4);
 }
 
 void obs_diagcorr::_start() {
@@ -76,68 +87,117 @@ void obs_diagcorr::_meas(bool loadedobs, const int& nthmeas) {
   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);
+	gsl_matrix_complex_set(measurements[nthmeas][itsep], icorr, jcorr,
+			       gsl_complex_rect( OM[icorr*4+jcorr][itsep].real(), OM[icorr*4+jcorr][itsep].imag() ));
       }
+  
+  for (int icorr=0; icorr<4; icorr++) {
+    gsl_vector_complex_set(measurements_disco[nthmeas], icorr,
+			   gsl_complex_rect( OM[16][icorr].real(), OM[16][icorr].imag() ) );
+  }
 };
 
+void obs_diagcorr::invertcm (gsl_matrix_complex *invm, const gsl_matrix_complex *m)
+{
+  int s;
+  gsl_permutation * perm = gsl_permutation_alloc (m->size1);
+  
+  gsl_matrix_complex *mtmp = gsl_matrix_complex_alloc(m->size1,m->size2);
+  gsl_matrix_complex_memcpy(mtmp, m);
+
+  gsl_linalg_complex_LU_decomp (mtmp, perm, &s);
+  gsl_linalg_complex_LU_invert (mtmp, perm, invm);
+
+  gsl_matrix_complex_free(mtmp);
+}
+
 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 *evec = gsl_matrix_complex_alloc(4,4);  
+  gsl_eigen_hermv_workspace *wspace = gsl_eigen_hermv_alloc(100);
+
+  gsl_eigen_hermv(m, v, evec, wspace);
+  gsl_eigen_hermv_sort(v, evec, GSL_EIGEN_SORT_VAL_ASC);
+  
+  gsl_eigen_hermv_free(wspace);
+  gsl_matrix_complex_free(evec);
 }
 
 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_matrix_complex *tmpmatrix2 = gsl_matrix_complex_alloc(4,4);
   gsl_vector *tmpvec = gsl_vector_alloc(4);
+  gsl_vector_complex *tmpvecc = gsl_vector_complex_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_complex *totaldisco = gsl_vector_complex_alloc(4);
+  gsl_matrix_complex *invC0 = NULL;
+  
+  gsl_vector_complex_set_zero(totaldisco);
+  for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
+    gsl_vector_complex_add( totaldisco, measurements_disco[imeas] );
+  
   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
+    gsl_matrix_complex_scale (tmpmatrix, gsl_complex_rect(1.0/O815->comargs.nmeas, 0.0) );
+    for (int icorr=0; icorr<4; icorr++)
+      for (int jcorr=0; jcorr<4; jcorr++) {
+	gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(totaldisco, icorr),
+						 gsl_complex_conjugate( gsl_vector_complex_get(totaldisco, jcorr) )
+						 );
+	discopart = gsl_complex_mul_real(discopart, pow(1.0/O815->comargs.nmeas,2));
+	gsl_matrix_complex_set( tmpmatrix, icorr, jcorr, 
+				gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ),
+						 discopart
+						 )
+				);
+      }
+
+    if (invC0 == NULL) {
+      invC0 = gsl_matrix_complex_alloc(4,4);
+      invertcm(invC0, tmpmatrix); 
+    }
+
+    gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
+		    invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
+    cdiag(jackres, tmpmatrix2);
 
     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 );
-
-      cdiag(tmpvec, tmpmatrix);
-
+      gsl_matrix_complex_scale ( tmpmatrix, gsl_complex_rect(1.0/(O815->comargs.nmeas-1), 0.0) );
+      gsl_vector_complex_memcpy( tmpvecc, totaldisco );
+      gsl_vector_complex_sub( tmpvecc, measurements_disco[imeas] );
+      for (int icorr=0; icorr<4; icorr++)
+	for (int jcorr=0; jcorr<4; jcorr++) {
+	  gsl_complex discopart = gsl_complex_mul( gsl_vector_complex_get(tmpvecc, icorr),
+						   gsl_complex_conjugate( gsl_vector_complex_get(tmpvecc, jcorr) )
+						   );
+	  discopart = gsl_complex_mul_real(discopart, pow(1.0/(O815->comargs.nmeas-1),2));
+	  gsl_matrix_complex_set( tmpmatrix, icorr, jcorr, 
+				  gsl_complex_sub( gsl_matrix_complex_get( tmpmatrix, icorr, jcorr ),
+						   discopart
+						   )
+				  );
+	}
+
+      gsl_blas_zgemm (CblasNoTrans, CblasNoTrans, gsl_complex_rect(1.0,0.0),
+		      invC0, tmpmatrix, gsl_complex_rect(0.0,0.0), tmpmatrix2);
+      
+      cdiag(tmpvec, tmpmatrix2);
       gsl_vector_sub(tmpvec, jackres);
       gsl_vector_memcpy(tmpvec2, tmpvec);
-      gsl_vector_mul(tmpvec, tmpvec2); //!!!!!!
-
+      gsl_vector_mul(tmpvec, tmpvec2);
       gsl_vector_add(jackerrornorm, tmpvec);
-    }
-    
+    } 
     gsl_vector_scale( jackerrornorm,
 		      (double)(O815->comargs.nmeas-1)/O815->comargs.nmeas );
 
@@ -148,24 +208,28 @@ void obs_diagcorr::_finish() {
 	   << "\t" << sqrt(gsl_vector_get(jackerrornorm,iev));
     }
     *out << endl;
-    
-    //for (int imeas=0; imeas<O815->comargs.nmeas; imeas++)
-      //gsl_matrix_
-      //  gsl_matrix_complex_memcpy( manymeans[imeas], precalc );
   }
 
   gsl_matrix_complex_free(totalval);
+  gsl_matrix_complex_free(tmpmatrix);
+  gsl_vector_free(tmpvec);
+  gsl_vector_complex_free(tmpvecc);
+  gsl_vector_free(tmpvec2);
+  gsl_vector_free(jackres);
+  gsl_vector_free(jackerrornorm);
+  gsl_vector_complex_free(totaldisco);
 };
 
 void obs_diagcorr::corrCompute()
 {
   complex<double> phislice[4][O815->comargs.lsize[1]];
 
+  for (int icorr=0; icorr<4; icorr++)
+    OM[16][icorr] = 0;
+  
   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 icorr=0; icorr<4; icorr++)
+      phislice[icorr][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 ];
@@ -173,13 +237,14 @@ void obs_diagcorr::corrCompute()
       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++) {
+      phislice[icorr][it] /= spatialV;
+      OM[16][icorr] += phislice[icorr][it];
+    }
   }
 
-  for (int icorr = 0; icorr < 4; icorr++)
+  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;
@@ -190,6 +255,9 @@ void obs_diagcorr::corrCompute()
 	OM[icorr*4+jcorr][itsep] /= O815->comargs.lsize[1];
       }
     }
+    
+    OM[16][icorr] /= O815->comargs.lsize[1];
+  }
 }
 
 #endif