u1casc-ordinary/obs_corr.hpp

   1 #ifndef OBS_CORR_HPP
   2 #define OBS_CORR_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 #define U1 0
  16 #define U2 1
  17 #define F1 2
  18 #define F2 3
  19
  20 using namespace std;
  21
  22 class obs_corr : public o815::obs {
  23
  24 public:
  25   obs_corr(o815 *_O815);
  26
  27 private:
  28   void _start();
  29   void _meas(bool loadedobs, const int& nthmeas);
  30   void _finish();
  31
  32   void corrCompute();
  33   static complex<double> effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para);
  34   static void preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para);
  35
  36   sim *Sim;
  37   obstat< complex<double>, complex<double>  > oC[4];
  38   obstat< complex<double>, complex<double>  > oD[4];
  39
  40   int spatialV;
  41
  42   complex<double> *OM[4];
  43 };
  44
  45 obs_corr::obs_corr(o815 *_O815) : o815::obs("corr",
  46                                             _O815->paraQ->getParaNames() +
  47                                             "tsep"
  48                                             ":u1_real:u1_imag:u1_abs:u1_effmass:u1_effmass_err"
  49                                             ":u2_real:u2_imag:u2_abs:u2_effmass:u2_effmass_err"
  50                                             ":f1_real:f1_imag:f1_abs:f1_effmass:f1_effmass_err"
  51                                             ":f2_real:f2_imag:f2_abs:f2_effmass:f2_effmass_err",
  52                                             _O815, 4*sizeof(complex<double>)*(_O815->comargs.lsize[1]/2+1) ) {
  53   for (int ivar = 0; ivar<4; ivar++)
  54     OM[ivar]  = (complex<double>*)( obsMem + ivar*sizeof(complex<double>)*( _O815->comargs.lsize[1]/2 + 1 ) );
  55
  56   Sim = (sim*)O815->Sim;
  57   spatialV = O815->comargs.lsize[0] * O815->comargs.lsize[0] * O815->comargs.lsize[0];
  58 }
  59
  60 void obs_corr::_start() {
  61   //*out << "OBS_test: start" << endl;
  62 };
  63
  64 void obs_corr::_meas(bool loadedobs, const int& nthmeas) {
  65   if (!loadedobs)
  66     corrCompute();
  67
  68   oC[U1].addMeas( OM[U1], O815->comargs.lsize[1]/2+1 );
  69   oC[U2].addMeas( OM[U2], O815->comargs.lsize[1]/2+1 );
  70   oC[F1].addMeas( OM[F1], O815->comargs.lsize[1]/2+1 );
  71   oC[F2].addMeas( OM[F2], O815->comargs.lsize[1]/2+1 );
  72 };
  73
  74 void obs_corr::_finish() {
  75   int compid_corr[O815->comargs.lsize[1]/2][4];
  76   int compid_effmass[O815->comargs.lsize[1]/2-1][4];
  77
  78   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++)
  79     for (int icorr = 0; icorr < 4; icorr++)
  80       compid_corr[itsep][icorr] = oC[icorr].computeEasy(itsep);
  81
  82
  83   for (int icorr = 0; icorr < 4; icorr++) {
  84     for (int itsep = 0; itsep < O815->comargs.lsize[1]/2-1; itsep++) {
  85       pair<int,int> effmasspass( itsep, O815->comargs.lsize[1]/2 );
  86       compid_effmass[itsep][icorr] = oC[icorr].computeJack(obs_corr::preEffMass, obs_corr::effMass, &effmasspass);
  87     }
  88   }
  89
  90   for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
  91     *out << O815->paraQ->getParaVals();
  92     *out << "\t" << itsep;
  93     for (int icorr = 0; icorr < 4; icorr++) {
  94       *out << "\t" << real( oC[icorr].getMean(compid_corr[itsep][icorr]) )
  95            << "\t" << imag( oC[icorr].getMean(compid_corr[itsep][icorr]) )
  96            << "\t" << abs ( oC[icorr].getMean(compid_corr[itsep][icorr]) );
  97
  98
  99       if ( itsep < O815->comargs.lsize[1]/2-1 )
 100         *out << "\t" << real( oC[icorr].getMean(compid_effmass[itsep][icorr]) )
 101              << "\t" << real( oC[icorr].getErr (compid_effmass[itsep][icorr]) );
 102       else
 103         *out << "\t" << NAN << "\t" << NAN;
 104     }
 105     *out << endl;
 106   }
 107
 108   for (int icorr = 0; icorr < 4; icorr++) {
 109     oC[icorr].reset();
 110     oD[icorr].reset();
 111   }
 112 };
 113
 114 void obs_corr::corrCompute()
 115 {
 116   complex<double> phislice[4][O815->comargs.lsize[1]];
 117
 118   OM[U1][O815->comargs.lsize[1]/2] = 0;
 119   OM[U2][O815->comargs.lsize[1]/2] = 0;
 120   OM[F1][O815->comargs.lsize[1]/2] = 0;
 121   OM[F2][O815->comargs.lsize[1]/2] = 0;
 122
 123   for (int it = 0; it < O815->comargs.lsize[1]; it++) {
 124     phislice[U1][it] = 0;
 125     phislice[U2][it] = 0;
 126     phislice[F1][it] = 0;
 127     phislice[F2][it] = 0;
 128
 129     for (int ix = 0; ix < spatialV; ix++) {
 130       phislice[U1][it] += conj(Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 131       phislice[U2][it] += conj(Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ]) * Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 132       phislice[F1][it] += Sim->phi[ 0*Sim->lsize4 + it*spatialV + ix ];
 133       phislice[F2][it] += Sim->phi[ 1*Sim->lsize4 + it*spatialV + ix ];
 134     }
 135     phislice[U1][it] /= spatialV;
 136     phislice[U2][it] /= spatialV;
 137     phislice[F1][it] /= spatialV;
 138     phislice[F2][it] /= spatialV;
 139
 140     OM[U1][O815->comargs.lsize[1]/2] += phislice[U1][it];
 141     OM[U2][O815->comargs.lsize[1]/2] += phislice[U2][it];
 142     OM[F1][O815->comargs.lsize[1]/2] += phislice[F1][it];
 143     OM[F2][O815->comargs.lsize[1]/2] += phislice[F2][it];
 144   }
 145
 146   for (int icorr = 0; icorr < 4; icorr++) {
 147     for (int itsep = 0; itsep < O815->comargs.lsize[1]/2; itsep++) {
 148       OM[icorr][itsep] = 0;
 149
 150       for (int it = 0; it < O815->comargs.lsize[1]; it++)
 151         OM[icorr][itsep] += phislice[icorr][ (it+itsep)%O815->comargs.lsize[1] ] * conj( phislice[icorr][it] );
 152
 153       OM[icorr][itsep] /= O815->comargs.lsize[1];
 154     }
 155
 156     OM[icorr][O815->comargs.lsize[1]/2] /= O815->comargs.lsize[1];
 157   }
 158 }
 159
 160 void obs_corr::preEffMass(vector< vector < complex<double> > > *allVals, vector < complex<double> > *preCalculated, void *para) {
 161   pair<int,int> *myparas = (pair<int,int>*)para;
 162
 163   preCalculated->push_back(0);
 164   preCalculated->push_back(0);
 165   preCalculated->push_back(0);
 166
 167   for(vector< vector< complex<double> > >::iterator valIt = allVals->begin(); valIt != allVals->end(); ++valIt) {
 168     (*preCalculated)[0] += (*valIt)[myparas->first];
 169     (*preCalculated)[1] += (*valIt)[myparas->first+1];
 170     (*preCalculated)[2] += (*valIt)[myparas->second];
 171   }
 172 }
 173
 174 complex<double> obs_corr::effMass(vector < complex<double> > *preCalculated, vector< complex<double> > *excludedMeas, int nmeas, void *para) {
 175   pair<int,int> *myparas = (pair<int,int>*)para;
 176
 177   double disconnected = norm( ( (*preCalculated)[2] - (*excludedMeas)[myparas->second] ) / (complex<double>)(nmeas-1) );
 178
 179   return std::log( abs(
 180                        ( ( (*preCalculated)[0] - (*excludedMeas)[myparas->first] ) / (complex<double>)(nmeas-1) - disconnected ) /
 181                        ( ( (*preCalculated)[1] - (*excludedMeas)[myparas->first+1] ) / (complex<double>)(nmeas-1) - disconnected )
 182                         ) );
 183 }
 184
 185 #endif