Added missing isospin correlators.

diff --git a/u1casc-ordinary/obs_corriso00.hpp b/u1casc-ordinary/obs_corriso00.hpp
new file mode 100644 (file)
index 0000000..ca4f2db
--- /dev/null
+++ b/u1casc-ordinary/obs_corriso00.hpp
@@ -0,0 +1,150 @@
+#ifndef OBS_CORRISO00_HPP
+#define OBS_CORRISO00_HPP
+
+#include "latlib/o815/o815.h"
+
+#include "latlib/writeout.h"
+
+#include "latlib/obstat.hpp"
+
+#include <iostream>
+#include <complex>
+
+#include <cmath>
+
+using namespace std;
+
+class obs_corriso00 : public o815::obs {
+
+public:
+  obs_corriso00(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;
+
+  int spatialV;
+
+  complex<double> *OM;
+};
+
+obs_corriso00::obs_corriso00(o815 *_O815) : o815::obs("corriso00", 
+                                                     _O815->paraQ->getParaNames() + 
+                                                     "tsep"
+                                                     ":iso00_real:iso00_imag:iso00_abs:iso00_mass:iso00_mass_err", 
+                                                     _O815, sizeof(complex<double>)*(_O815->comargs.lsize[1]/2+1) ) {
+  
+  OM  = (complex<double>*)(obsMem);
+  
+  Sim = (sim*)O815->Sim;
+  spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
+}
+
+void obs_corriso00::_start() {
+  //*out << "OBS_test: start" << endl;
+};
+
+void obs_corriso00::_meas(bool loadedobs, const int& nthmeas) {
+  if (!loadedobs)
+    corrCompute();
+  
+  oC.addMeas( OM, O815->comargs.lsize[1]/2+1 );
+};
+
+void obs_corriso00::_finish() {
+  int compid_corr[O815->comargs.lsize[1]/2];
+  int compid_effmass[O815->comargs.lsize[1]/2-1];
+
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++)
+    compid_corr[itsep] = oC.computeEasy(itsep);
+  
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) {
+    pair<int,int> effmasspass( itsep, O815->comargs.lsize[1]/2 );
+    compid_effmass[itsep] = oC.computeJack(obs_corriso00::preEffMass, obs_corriso00::effMass, &effmasspass);
+  }
+  
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+    *out << O815->paraQ->getParaVals();
+    *out << "\t" << itsep;
+   
+    *out << "\t" << real( oC.getMean(compid_corr[itsep]) )
+        << "\t" << imag( oC.getMean(compid_corr[itsep]) )
+        << "\t" << abs ( oC.getMean(compid_corr[itsep]) );
+      
+    if ( itsep < O815->comargs.lsize[1]/2-1 )
+      *out << "\t" << real( oC.getMean(compid_effmass[itsep]) )
+          << "\t" << real( oC.getErr (compid_effmass[itsep]) );
+    else
+      *out << "\t" << NAN << "\t" << NAN;
+    
+    *out << endl;
+  }
+  
+  oC.reset();
+}
+
+
+void obs_corriso00::corrCompute()
+{
+  complex<double> phislice[O815->comargs.lsize[1]];
+  
+  OM[O815->comargs.lsize[1]/2] = 0;
+
+  for (int it = 0; it < O815->comargs.lsize[1]; it++) {
+    phislice[it] = 0;
+    
+    for (int ix = 0; ix < spatialV; ix++)
+      phislice[it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]
+       + conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
+      
+    phislice[it] /= sqrt(2)*spatialV;
+
+    OM[O815->comargs.lsize[1]/2] += phislice[it];
+  }
+
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+    OM[itsep] = 0;
+
+    for (int it = 0; it < O815->comargs.lsize[1]; it++)
+      OM[itsep] += phislice[ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[it] );
+      
+    OM[itsep] /= O815->comargs.lsize[1];
+  }
+
+  OM[O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1];
+}
+
+void obs_corriso00::preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para) {
+  pair<int,int> *myparas = (pair<int,int>*)para;
+  
+  preCalculated->push_back(0);
+  preCalculated->push_back(0);
+  preCalculated->push_back(0);
+
+  for(vector< vector< complex<double> > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
+    (*preCalculated)[0] += (*valIt)[myparas->first];
+    (*preCalculated)[1] += (*valIt)[myparas->first+1];
+    (*preCalculated)[2] += (*valIt)[myparas->second];
+  }
+}
+
+complex<double> obs_corriso00::effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para) {
+  pair<int,int> *myparas = (pair<int,int>*)para;
+
+  double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex<double>)(nmeas-1) );
+  
+  return std::log( abs( 
+                      ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex<double>)(nmeas-1) - disconnected ) / 
+                      ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex<double>)(nmeas-1) - disconnected ) 
+                       ) );
+}
+
+#endif

diff --git a/u1casc-ordinary/obs_corriso10.hpp b/u1casc-ordinary/obs_corriso10.hpp
new file mode 100644 (file)
index 0000000..3bdb7b5
--- /dev/null
+++ b/u1casc-ordinary/obs_corriso10.hpp
@@ -0,0 +1,150 @@
+#ifndef OBS_CORRISO10_HPP
+#define OBS_CORRISO10_HPP
+
+#include "latlib/o815/o815.h"
+
+#include "latlib/writeout.h"
+
+#include "latlib/obstat.hpp"
+
+#include <iostream>
+#include <complex>
+
+#include <cmath>
+
+using namespace std;
+
+class obs_corriso10 : public o815::obs {
+
+public:
+  obs_corriso10(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;
+
+  int spatialV;
+
+  complex<double> *OM;
+};
+
+obs_corriso10::obs_corriso10(o815 *_O815) : o815::obs("corriso10", 
+                                                     _O815->paraQ->getParaNames() + 
+                                                     "tsep"
+                                                     ":iso10_real:iso10_imag:iso10_abs:iso10_mass:iso10_mass_err", 
+                                                     _O815, sizeof(complex<double>)*(_O815->comargs.lsize[1]/2+1) ) {
+  
+  OM  = (complex<double>*)(obsMem);
+  
+  Sim = (sim*)O815->Sim;
+  spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
+}
+
+void obs_corriso10::_start() {
+  //*out << "OBS_test: start" << endl;
+};
+
+void obs_corriso10::_meas(bool loadedobs, const int& nthmeas) {
+  if (!loadedobs)
+    corrCompute();
+  
+  oC.addMeas( OM, O815->comargs.lsize[1]/2+1 );
+};
+
+void obs_corriso10::_finish() {
+  int compid_corr[O815->comargs.lsize[1]/2];
+  int compid_effmass[O815->comargs.lsize[1]/2-1];
+
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++)
+    compid_corr[itsep] = oC.computeEasy(itsep);
+  
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) {
+    pair<int,int> effmasspass( itsep, O815->comargs.lsize[1]/2 );
+    compid_effmass[itsep] = oC.computeJack(obs_corriso10::preEffMass, obs_corriso10::effMass, &effmasspass);
+  }
+  
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+    *out << O815->paraQ->getParaVals();
+    *out << "\t" << itsep;
+   
+    *out << "\t" << real( oC.getMean(compid_corr[itsep]) )
+        << "\t" << imag( oC.getMean(compid_corr[itsep]) )
+        << "\t" << abs ( oC.getMean(compid_corr[itsep]) );
+      
+    if ( itsep < O815->comargs.lsize[1]/2-1 )
+      *out << "\t" << real( oC.getMean(compid_effmass[itsep]) )
+          << "\t" << real( oC.getErr (compid_effmass[itsep]) );
+    else
+      *out << "\t" << NAN << "\t" << NAN;
+    
+    *out << endl;
+  }
+  
+  oC.reset();
+}
+
+
+void obs_corriso10::corrCompute()
+{
+  complex<double> phislice[O815->comargs.lsize[1]];
+  
+  OM[O815->comargs.lsize[1]/2] = 0;
+
+  for (int it = 0; it < O815->comargs.lsize[1]; it++) {
+    phislice[it] = 0;
+    
+    for (int ix = 0; ix < spatialV; ix++)
+      phislice[it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]
+       - conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
+      
+    phislice[it] /= sqrt(2)*spatialV;
+
+    OM[O815->comargs.lsize[1]/2] += phislice[it];
+  }
+
+  for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
+    OM[itsep] = 0;
+
+    for (int it = 0; it < O815->comargs.lsize[1]; it++)
+      OM[itsep] += phislice[ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[it] );
+      
+    OM[itsep] /= O815->comargs.lsize[1];
+  }
+
+  OM[O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1];
+}
+
+void obs_corriso10::preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para) {
+  pair<int,int> *myparas = (pair<int,int>*)para;
+  
+  preCalculated->push_back(0);
+  preCalculated->push_back(0);
+  preCalculated->push_back(0);
+
+  for(vector< vector< complex<double> > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
+    (*preCalculated)[0] += (*valIt)[myparas->first];
+    (*preCalculated)[1] += (*valIt)[myparas->first+1];
+    (*preCalculated)[2] += (*valIt)[myparas->second];
+  }
+}
+
+complex<double> obs_corriso10::effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para) {
+  pair<int,int> *myparas = (pair<int,int>*)para;
+
+  double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex<double>)(nmeas-1) );
+  
+  return std::log( abs( 
+                      ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex<double>)(nmeas-1) - disconnected ) / 
+                      ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex<double>)(nmeas-1) - disconnected ) 
+                       ) );
+}
+
+#endif

diff --git a/u1casc-ordinary/u1casc-ordinary.cpp b/u1casc-ordinary/u1casc-ordinary.cpp
index 4803fe11ef74834e6116abcffb071edf7da033d4..b6b69c077cf01d694d81845b601ee962409e3f13 100644 (file)
--- a/u1casc-ordinary/u1casc-ordinary.cpp
+++ b/u1casc-ordinary/u1casc-ordinary.cpp
@@ -17,6 +17,8 @@ sim *Sim;
 #include "obs_corrphichi.hpp"
 #include "obs_corrphi.hpp"
 #include "obs_corrchi.hpp"
+#include "obs_corriso10.hpp"
+#include "obs_corriso00.hpp"
 
 o815::comoption specOps[] = {
   { "kappaone", required_argument, NULL, 'r', "set inverse mass kappa_1", "min:max:inc" },
@@ -95,6 +97,14 @@ void parseLonelyArgs()
       *O815->out->log << "MASTER: registered observable: corrchi" << endl << flush;
       O815->observables.push_back(new obs_corrchi(O815));
     }
+    else if ( strcmp(*lonit, "corriso10") == 0 ) {
+      *O815->out->log << "MASTER: registered observable: corriso10" << endl << flush;
+      O815->observables.push_back(new obs_corriso10(O815));
+    }
+    else if ( strcmp(*lonit, "corriso00") == 0 ) {
+      *O815->out->log << "MASTER: registered observable: corriso00" << endl << flush;
+      O815->observables.push_back(new obs_corriso00(O815));
+    }
   }
 }